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Introduction 


The study of the absorption, emission, and scattering of electromag- 
netic radiation as it passes through a medium with which it interacts is 
a fascinating subject involving the close interconnection of many disci- 
plines in mathematics and physics. This subject originated in the study 
of the properties and fate of radiant energy as it traverses stellar inte- 
riors, and much of the terminology and many of the definitions reflect 
the impetus given by the early researchers in the field. More recently, 
a great body of this theory has been applied to the study of the pas- 
sage of solar and terrestrial radiation through the Earth’s atmosphere, 
as well as to the study of radiation in the atmospheres of the other 
planets. In particular, studies of climate and climate models, and the 
detection and measurement of the distribution of water vapor, trace 
gases, and aerosols in the atmosphere have given additional importance 
to this topic, and literally hundreds of technical papers have been writ- 
ten in the past 20 or so years in which applications of radiative transfer 
(RT) theory have been made to these and other topics in atmospheric 
physics. 

The lack of standardization of symbols and terms in current radia- 
tive transfer literature has caused some difficulty, especially for the neo- 
phyte researcher, in comparing analyses and numerical results among 
the published texts and papers in radiative transfer theory; this conse- 
quently presents the new researcher with some difficulties in developing 
an integrated picture of, or feel for, this most fascinating subject. The 
present monograph is an attempt to alleviate this frustrating circum- 
stance by developing some of the fundamental concepts in RT theory, 
and by defining some of the more useful approximate solutions to the 
radiative transfer equation (RTE) using as consistent a set of definitions 
and symbols as is practical. This will hopefully make the newcomer’s 
transition to the more formal technical literature somewhat less painful. 

The radiative transfer equation appears in many forms in the 
literature, depending on the discipline, the area of application, and the 
whims of the writer. The various forms derived herein are those most 
generally encountered in atmospheric applications. As in any scientific 
discipline, the technical literature is generally written by experts for 
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experts, and, as a consequence of long familiarity with the basic theory, 
a great deal is generally omitted from their papers as being well known 
or implied, causing still more confusion to the researcher new to the 
field. Frequently, for example, one paper presents specialized forms 
of the RTE which supposedly represent the same physical situation 
as in another paper, and yet the physical forms of the corresponding 
equations are dissimilar. The present book will hopefully aid the reader 
in recognizing these differences and the reasons for them, and thus allow 
the reader to construct a mental link between the seemingly different 
results. 

Some of the classical solutions to the various forms of the RTE 
will be derived in detail. These will include the thin-atmosphere 
approximation, the single-scattering solution, various forms of the two- 
stream solutions, the Eddington solution, and the discrete ordinates 
method of Chandrasekhar. In some cases, numerical examples will be 
given so that the reader can develop a feel for the order of magnitude of 
the numbers involved. Appropriate caveats will be rendered concerning 
regions of applicability of the approximation methods. 

A working group of the Radiation Commission of the International 
Association of Meteorology and Atmospheric Physics (lAMAP) headed 
by Jacqueline Lenoble of the University of Lille, France, has edited 
an extremely comprehensive but very compact two-volume set of notes 
containing descriptions of all the presently used methods for comput- 
ing the radiative transfer through scattering atmospheres. Because of 
its scope, this document is difficult to use as a tutorial guide, but is 
an excellent reference source for the experienced researcher. Its title 
is “Standard Procedures to Compute Atmospheric Radiative Trans- 
fer in a Scattering Atmosphere”; it is published by the lAMAP, and 
is obtainable from NCAR, Boulder, Colorado. This document dis- 
cusses all the current problems in radiative theory, all the methods 
currently in fashion, and gives hundreds of references. It is highly rec- 
ommended for source material once the fundamentals of the present 
text are fully grasped. Prof. Lenoble hcts edited and revised a set of 
these documents, which is available as Radiative Transfer in Scatter- 
ing and Absorbing Atmospheres: Standard Computational Procedures^ 
A. Deepak Publishing, Hampton, Virginia, 1985 (ISBN 0-937194-05-0). 

The present book is not intended to be a textbook on radiative 
transfer theory, nor is it intended to be authoritative or complete — 
the author has neither the inclination nor the expertise to attempt 
such a monumental task. It is meant rather to be a set of notes, 
mathematically more detailed than one usually finds in a textbook 
or formal paper, presenting the derivations and solutions to various 
forms of the integro-differential equation which describes the transfer of 
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radiant energy through an absorbing, scattering, and emitting medium. 
The basic thrust of these notes is twofold: first, to provide the reader 
with a firm physical foundation of the basics of radiative transfer which 
will permit a ready transition to the more formal literature from which 
this foundation can be expanded; second, to present some of the more 
elementary, but perhaps more useful, solutions to the RTE in sufficient 
detail for the reader to be able to concentrate on the physical principles 
involved in these developments rather than being bogged down by a lot 
of superfluous mathematical detail, thereby helping the reader develop 
a physical feel for the way the various components interact and for their 
relative importance. 

To give proper credit where it is due, it should be mentioned that 
except for the speciflc papers referenced in the body of the text, most of 
the material for this monograph was extracted from three texts: those 
of Chandrasekhar (1960), Liou (1980), and Ozisik (1973) — in particular, 
chapters 1-5 of Chandrasekhar, chapters 1 and 6 of Liou, and chapters 
1, 8, and 9 of Ozisik. 
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Chapter 1 

Introductory Concepts 


All substances continuously emit electromagnetic radiation as a 
result of the thermal motions of the molecules and atoms of which 
they are made. The thermal agitation of these particles increases 
with temperature, and consequently, the emitted radiation frequencies 
increase with temperature. The wavelengths of these radiations range 
from several kilometers for very long wavelength radio transmissions 
down to 10”^^ cm and less for cosmic rays and beyond. A very rough 
and somewhat arbitrary division of the electromagnetic spectrum is 
given below in table 1-1. 

TABLE 1-1. THE ELECTROMAGNETIC SPECTRUM 


Wavelength 

Type of radiation 

o 

1 

o 

o 

o 

n 

B 

Radio, radar, TV, etc. 

10“^ to 10“^ cm 

Infrared 

10“^ to 10“^^ cm 

Visible 

10”® to 10”^ cm 

Ultraviolet 

to 10”® cm 

X-rays 

10“^^ to 10”^ cm 

Gamma rays 

? to 10“^^ cm 

Cosmic rays 


The term thermal radiation is normally reserved for radiation that 
can be detected as either heat or light, and so is generally applied to 
that region of the spectrum ranging from about 10“^ to 10“^ cm; i.e., 
the visible and infrared portions of the spectrum. Since we shall be 
primarily concerned with the infrared portion of the spectrum, a unit 
called the micron, equal to 10”"^ cm (or 10“^ m) will be used throughout 
these notes. In these terms, the thermal radiation regime extends from 
about 0.1 to 1000 microns. 

PRECtDiNG PAGE BLA.NK NOT RLWa^ 5 
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Now, in addition to emitting their own radiation, most atmospheric 
constituents, from molecules to water vapor droplets and aerosols, also 
affect radiation incident on them by the processes of absorption and 
scattering. Scattering can be thought of as the process of changing 
the direction of the incident radiation — in some cases, by changing 
the frequency of the scattered radiation. The problem of scattering 
radiation without changing its frequency is the main topic of these 
notes. 

The relative proportion of the incident radiation which is scattered 
in various directions by the scattering particles is a function of basically 
two parameters: the size of the particle relative to the wavelength of 
the incident radiation, and the (complex) index of refraetion of the 
material of which the particle is made. The shape of the particle is also 
quite important, but the theory developed thus far, the Mie theory, is 
reasonably complete only for spherical particles, although particles of 
cylindrical and flat plate geometries have recently been addressed with 
some success. The basic parameters resulting from these analyses are 
the particle’s phase function^ which describes the spatial distribution of 
the scattered radiation, the scatter cross section, which determines the 
fraction of the incident radiation which is scattered, and the absorption 
cross section, which, for particles with a nonzero imaginary index of 
refraction, defines the fraction of the incident radiation that is absorbed 
by the particle. 

The determination of these parameters is a subject of its own and 
will not be addressed in these notes — here, these parameters will be 
assumed to be known. The text by Liou, cited earlier, gives a good 
introduction to this subject, and the classical texts of van de Hulst 
(1957) and Deirmendjian (1969) should be consulted for more details. 
The texts by Kerker (1969) and Stratton (1941) are also quite readable 
and useful, and the excellent review paper by Hansen and Travis (1974) 
covers the scattering problem very concisely, as well as many of the 
other topics presented in the present text. 

The theory of radiant energy propagation can generally be consid- 
ered from two different viewpoints: classical electromagnetic wave the- 
ory, and quantum mechanics. 

The classical theory begins with Maxwell’s equations and considers 
the energy propagation characteristics of electromagnetic waves. How- 
ever, the classical theory generally ignores the microscopic interactions 
of the radiation with matter, and treats only the macroscopic behav- 
ior. Consequently, many of the parameters of interest in the study of 
the propagation of radiation through absorbing, scattering, and emit- 
ting media are defined quantities which must be determined through 
experiment. 
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The situation is quite similar to that found in classical thermody- 
namics. Acceptance of the first and second laws of thermodynamics, 
an introduction of the perfect gas law, and the concept of entropy al- 
low a great many mathematical statements to be made which correctly 
identify basic trends and gross features of thermodynamic systems in 
equilibrium. These classical concepts by themselves, however, generally 
do not permit the detailed calculation of numerical results. Certain con- 
cepts and groups of parameters are related to others through arbitrary 
constants of proportionality which must be experimentally determined. 
Such parameters as specific heats, heat transfer, diffusion coefficients, 
thermal conductivity, and viscosity coefficients are merely “constants 
of proportionality,” and classical thermodynamics offers no means of 
directly computing their numerical magnitudes from first principles or 
of predicting the way in which these parameters will vary with such 
macroscopic thermodynamic properties as temperature, pressure, etc. 

Classical statistical mechanics does attack these problems within the 
framework of classical physics by making some hypotheses concerning 
the molecular structure of the material; i.e., it assumes a mathematical 
“model” of the system. In this way, many of the above-mentioned 
coefficients can be computed in terms of the modeled physical properties 
of the molecules and the local properties of the system. These results, 
which are to a greater or lesser extent constrained by the fidelity 
of the assumed model, generally can predict the gross characteristics 
of these coefficients adequately and, when applied to systems which 
are known to fall within the realm of “classical physics,” predict 
numerical magnitudes reasonably well. However, ultimately, an appeal 
to quantum statistical mechanics must be made to account for behavior 
which classical theory cannot handle. Unfortunately, the mathematical 
structure of these equations is generally very complex, and much of the 
insight offered by the classical theory is lost. 

Initially, we shall adopt the quantum mechanical approach for the 
analysis of the radiation field. That is, we shall consider the field 
to be composed of photons rather than waves, and shall define the 
basic properties of the field in these terms. However, frequently 
an appeal to the classical approach will be made in the interest of 
clarity or expediency. For example, while the basic property of the 
radiation field, the spectral intensity, will be defined from the photon 
model, the concepts of absorption and emission coefficients will be 
introduced in classical terms as constants of proportionality in equations 
which describe the changes in spectral intensity as the radiation passes 
through and interacts with an optically active medium. The two 
approaches will also be combined in the derivation of the basic radiative 
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transfer equation, where the absorption and emission coefficients will 
be related to the annihilation and production of photons. 

The concepts of absorption and emission coefficients can be devel- 
oped by a formal quantum mechanical approach which relates the ab- 
sorption and emission properties of the medium to the Einstein Transi- 
tion Probability coefficients, which ultimately permit their calculation 
in terms of the microscopic properties of the medium. 

Definitions 

The initial concepts to be presented below are those which for the 
writer were initially the most difficult to understand and to relate to 
physically meaningful ideas. Consequently, extreme and frequently 
painful notational rigor will be adhered to initially, and parenthetical 
references to variable dependencies will be used in abundance. Both 
of these cumbersome nuisances will be relaxed or dropped in subse- 
quent sections, as familiarity with their concepts and implications will 
hopefully have been attained by then. 

Suppose one has an arbitrary volume, which contains N photons. 
These photons all travel with the speed of light, c, but they have 
definite distributions of energy and directions of motion. If the volume 
is ^LSSumed to be in thermodynamic equilibrium, the energy distribution 
is given by the well-known Planck function, 

B m 

' c2[exp(/ij//A:r) - 1] 

in which h is Planck^s constant, 6.626 x joule-sec, c is the velocity 

of light, 2.998 x 10® m/sec, k is Boltzmann's constant, 1.381 x 10”^® 
joule/deg, ly is the frequency in hertz, and T is the absolute temperature 
in kelvins. Then, Bi^{T) has dimensions of watts/(m^-sec-st). We also 
follow Planck in postulating that if a particular photon has a frequency 
associated with it, then the energy of the photon is hu. 

Now consider the quantity n, where n — N/V is the total number of 
photons per unit volume, with all permissible energies and traveling in 
all directions. Of all these photons, let us single out all the ones whose 
energies lie in the range hu to h{iy -{- di/), and let symbolize these 
selected photons. Obviously, then, 

r<x> 

n = I Hu dv (1-1) 

Let us further restrict our selection of photons to all of the rij^ photons 
which are traveling in a specified direction defined by the unit vector n, 


8 


Chapter 1 


and which lie in a differential solid angle centered on fl. (See fig. 1-1.) 
From this, we define a photon distribution function^ fuy as the number of 
photons per unit volume having the direction of propagation fl within 
the solid angle dU, whose energies lie in the range of hu to h{u + du), 
and which are passing through a unit area in a unit time. Then, 

= y U{Cl) dn (1-2) 



Figure 1-1. Solid angle and direction of travel of the selected photons. 

Now, if we consider the area element dA whose normal n makes 
an angle 9 with the fl-vector, then dA cos 6 is the projected area of 
dA normal to the direction of propagation n. If the photons are 
traveling with a velocity c then in time dt the total volume enclosing 
all the photons which have passed through dA in the direction of fl 
is {dA cos 9) (c dt), and thus, the number of selected photons in this 
volume is cf {Cl) cos 9 dA dt du dCt. Since each photon has an energy 
fti/, the total energy of all the selected photons is 

dWt, = hi^cU{n)cos9 dA dt du dCl (1-3) 

From this basic expression we can extract all of the definitions we need 
for our development, and hopefully for understanding other writers’ 
definitions. 

First of all 


= {dAcos'^dt d. dh = 

is defined as the spectral intensity or radiance^ and is, at least in 
theoretical developments, perhaps the most fundamental and useful 
property of a radiation field. It can be seen that the radiance is defined 
as the total energy per unit time in the frequency interval di/ crossing 
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the unit projected surface area normal to the direction of propagation 
and in the infinitesimal solid angle centered around the direction of 
propagation. 

Also from equation (1-3) we write 


dpi/ — 


dWi/ 

{dAcos0){c dt)du 


= hucfi/{Vi) df] 


and so 

Pu = J hi^uin) dn (1-5) 

is the total energy per unit volume of all the photons whose energy 
range is hi^ to h{iy + di/), but which are traveling in any direction. This 
is called the spectral energy density, and can be related to the radiance, 
by the use of equation (1-4) 

Pi. = ^ J dU ( 1 - 6 ) 

We now accept Ii/{ft) as our fundamental parameter, write equa- 
tion (1-3) in terms of this parameter, and use 

dWi/ = Ii/{fl) dAcosO dt dU du (1-7) 

as our basic equation. This equation is frequently presented as an 
intuitive relation, relating the total energy functionally to the area 
element, frequency, and the direction of propagation, in which case the 
radiance is frequently inserted as merely a constant of proportionality — 
hardly an auspicious introduction for such an important parameter. 
It can be seen from equation (1-4), however, that the radiance can 
be defined from more fundamental principles, and has a real physical 
identity of its own. 

The quantity 


= dA Tin (‘-s) 

is called the spectral emissive power. This is the total energy per unit 
time in the frequency interval du crossing the total unit area into the 
unit solid angle centered about the direction of propagation fl, and is 
a function of 6, as distinct from the definition of Ii/{Cl). 

From equations (1-7) and (1-8) the quantity 

^ dn 

dA dt du 
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or 


Fu 


j cos$ dQ 


(1-9) 


is called the spectral flux or irradiance^ and is probably the second most 
useful property describing the radiation field. This is the total energy 
within the frequency range du^ passing through the unit area per unit 
time, traveling in all possible directions. 

Lastly, we define 


dE = 


dWu 
dA dt 


= Ity{fl)cosO dQ du 


or 





cos 6 dU du = 



( 1 - 10 ) 


as the total emissive power. This is the total energy, or total flux per 
unit volume at all frequencies and in all directions passing through the 
unit area in unit time. 

This completes the set of basic definitions. It is hoped that by 
appealing to the corpuscular approach, rather than the classical concept 
of waves with their associated energies and intensities, the above 
definitions will be easier to grasp. 

As somewhat of an important aside, let us now make the assumption 
that all directions of motion of the photons are equally probable. 
This defines the concept of isotropy of radiation, and is one of the 
characteristics of blackbody radiation. All of the definitions given thus 
far are perfectly general and apply to any radiation field. In an isotropic 
field, some of these appear in a simpler algebraic form which will 
perhaps assist the reader to recall their physical significance. 

From equation (1-2), njy = or 


U = (HI) 

47T 

or, the number of photons per unit volume in the direction of any solid 
angle is equal to the total number of photons per unit volume divided 
by the area of the unit sphere; i.e., the directions of travel are uniformly 
distributed over the unit sphere. From equation (1-6) 


Pv = 


47r/t/ 

c 


and using equations (1-4) and (1-11) 


Pi, = 4nhufiy{ft) = huuiy 


( 1 - 12 ) 


(1-13) 
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The energy density of all photons in the frequency interval du traveling 
in any direction is equal to the energy per photon, times the number 


of selected photons. 



From equation (1-9) 



Fi/ = y' COS0 dU 

r2iv r7r/2 

= Ii/ 1 / COS 9 sin 9 d9 d(j) 

Jo Jo 

(1-14) 


II 

(1-15) 

and from equation (1-10) 

E = 7tI 

(1-16) 


where / = du is called the total intensity^ so that the total 

emissive power E bears the same relationship to the total intensity I as 
the radiant flux Fi, bears to the radiance Note in equations (1-4) and 
(1-14) that the upper limit in the 6 integral is tt/2 rather than tt. This is 
because these quantities are usually defined in the literature relative to 
an emitting surface^ and hence can only emit into a hemisphere centered 
on the elemental area. Thus Fi^ is sometimes called the hemispherical 
spectral radiant flux— quite a mouthful— and E is sometimes called the 
hemispherical total emissive power. Note, however, that the net flux is 
found by integrating over the whole unit sphere. 

Substituting equation (1-12) into equation (1-16) 

E^^f p^du (1-17) 

and using equation (1-13) 

E = ^ J huuu du (1-18) 

Finally, using equation (1-12) to eliminate I in equation (1-15) 

Fjy ~ (1-19) 

which relates the spectral radiant flux to the spectral energy density. 

Up to this point, we have used a somewhat quantum mechanical 
approach, in that we have considered the radiation field to be composed 
of photons rather than waves, as in the classical approach. All of the 
previous equations and their relationships to one another could have 
been derived from classical electromagnetic theory, naturally, and in 
fact historically have been derived in just this manner. 
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We will now shift temporarily to a more macroscopic approach in 
the sense that we will define some quantities which are used to describe 
the way radiation energy interacts with the medium through which it 
is propagating in terms of experimentally derived characteristics rather 
than formal analysis. We will find that certain observed properties of 
the radiation field are related to certain properties of the medium, and 
prescribe certain constants of proportionality to form useful mathemat- 
ical relationships. Classically, these constants are found by experiment, 
much the same way some of the classical thermodynamic coefficients 
mentioned earlier are found. Quantum mechanically, they can, at least 
in theory, be derived from considerations of the molecular structure 
of the medium and the electromagnetic interactions between the field 
produced by the medium and that produced by the radiation. 

Absorption of Radiation 

Consider a beam of monochromatic radiation of specific intensity 
/^y(r, n'); note, we now include the spatial dependence r — confined to 
an element of solid angle dXV that is incident normal to the surface 
dA of a slab of optically active material of thickness ds. (See fig. 1- 
2.) As the radiation passes through the slab, some of the photons 
will be absorbed by the material in the slab, some will be scattered 
out of the beam by the material, and the rest will emerge from the 
opposite face of the slab. We confine ourselves for the present to the 
photons which are absorbed. Define Ki^{r) as the spectral volumetric 
absorption coefficient. This coefficient has units of and represents 
the fraction of the incident radiation that is absorbed by the matter in 
the slab per unit length along the path of the incident radiation. Then, 
the total amount of radiation absorbed by the slab per unit time, per 
unit frequency interval in the solid angle dU^ is 

K^{r)I^{r, n') dn' dA ds (1-20) 


il') 



Figure 1-2. Infinitesimal cone of radiation intensity impinging on a thin slab. 

It can be shown (e.g., Sparrow and Cess, pp. 17-18) that l/Ku{r) 
can be interpreted as the mean free path for photon absorption; i.e., 1/e 
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of the incident photons will be absorbed within a distance of 1/Kiy{r) 
of the front surface. 

The volumetric absorption coefficient can be related to the more 
commonly used molecular absorption coefficient^ and the mass 

absorption coefficient, K^{v), as follows: assume that the optically 
active material in the slab has a number density of nm{^) molecules/m^. 
Each molecule has an absorption cross section of m^/molecule 

associated with it. Then the total absorption of the incident radiation 
in the length ds will be 

X^(r)/^(r, n') dn' n^(r) dA ds (1-21) 

Comparison of equation (1-20) with equation (1-21) reveals that 

K^{v) = (1-22) 

Similarly, if the optically active material has a mass density of 
Pm hg/m^, then the analogue to equation (1-21) is 

K^{t)Iu{t, n') dU' pm{T) dA ds (1-23) 

and 

K,{r) = Kt{r)pm{r) (1-24) 

Note that the units for KJP are m^/molecule, and for are m^/kg. 

Scattering of Radiation 

In addition to the attenuation of the incident beam by absorption, 
some of the photons of the incident beam are removed by the process 
of scattering. Let denote the spectral volumetric scattering 

coefficient. This coefficient has dimensions of m”^, and represents the 
fraction of the incident radiation that is scattered by the optically active 
material in the slab, in all directions, per unit length in the slab. Thus, 
the quantity 

au{r)I^{r,n') dn' (1-25) 

is the amount of incident radiation scattered in a unit length by the 
matter in all directions, per unit time and per unit frequency centered 
about 2 /. 

This relation does not supply any information about the directional 
distribution of the scattered radiation. We therefore introduce the 
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concept of the phase function, such that 

^p,(n,n')dn (1-26) 

Air 

describes the probability that the incident radiation, 7i/(r, n'), will be 
scattered from the solid angle dO' centered about fl' into an element of 
solid angle dQ centered about the direction of Q. The factor An is the 
total solid angle, and is introduced for normalization 

^ /p^(n,n')dn = i (1-27) 

Air J 

n 

which says that all of the scattered radiation must go somewhere in the 
unit sphere. 

We should note that some authorities, notably Chandrasekhar, 
define the integral in equation (1-27) as 


1 

Air 


j p^(n,n') dn = cj^ 

n 


where ojiy is the single-scattering albedo^ a concept to be introduced 
later, and thus represents the fraction of the total incident energy lost 
from the beam due to scattering only. Many authorities currently follow 
this practice; nonetheless, more and more experts seem to be adopting 
equation (1-27), which is preferable, in this writer’s opinion, as it allows 
the parameter to be injected into the formal Radiative Transfer 
Equation (RTE) somewhat less artificially. Hence equation (1-27) is 
the definition used in the next chapter when deriving the RTE. 

Putting equations (1-25) and (1-26) together, then, 


[a^(r)/^(r,n') dfl') 


Air 


p^(n,n') dn 


is the amount of the incident radiation which is scattered by the slab 
per unit time, volume, etc., into an element of solid angle dfl centered 
about n. Integrating this expression over all angles of incidence gives 

dn f iu{r,n') Pu{n,ii') dn' (1-28) 

Air J 

O' 
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which is the total radiation scattered into the element dfi from all 
directions of incidence per unit time, etc. 

Figure F3 shows the geometry of the scattering process. The angle 
6q is called the scattering angle (what else?). One usually assumes that 
the phase function depends only on the scattering angle. In that case, 
one usually writes 

p^{n,n') = p^{cos0o) (1-29) 

and hence writes equation (1-28) as 
1 f2n f7r/2 

— <7i^(r) dfl / /iy(r, 0',0')P^(cos0o) d<t/ (1-30) 
47t Jo ^0 

where the angles 9^ and are the usual colatitude and azimuthal polar 
coordinates. (See fig. 1-4.) Letting 9 and (j) represent the corresponding 
quantities for the scattered ray, we can write the expression for the unit 
vectors ft and ft\ 



sin 9 cos (f) 


sin 9^ cos (p' 

n = 

sin 0 sin 0 

n' = 

sin 9^ sin (p^ 


cos 9 


cos 9^ 


Then, since cos0q ~ O ♦ H', 

cos 9q — cos 9 cos 9^ 4- sin 9 sin 9^ cos{4>^ — (j>) 


(1-31) 


(1-32) 
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Figure 1-4. Spherical coordinate system used to define the scattering angle. The 

Z-3X1S is normal to the slab. 

or, as it is usually written, with /i = cos 9 

fiQ = ^ ^ cos(0' - 0) (1-33) 

It might be worthwhile to elaborate somewhat on the form of the 
phase function as given by equation (1-29), even though some of the 
ideas we use will not be formally introduced until a later section. We are 
interested here in the form of the phase function when it is a function 
of the scattering angle only, and not a function of azimuth. This is a 
constraint that is almost universally applied in the literature. 

For most atmospheric applications, the phase function has a shape 
which generally resembles the sketch in figure 1-5 — the figure is rota- 
tion ally symmetric. The scatter function generally has a small backscat- 
ter component (a), one or more “side-lobes” of various angular arrange- 
ments and magnitudes (b), and generally a strongly forward-scattering 
peak (c). The ratio of forward to backward scattering may in many 
cases exceed several hundreds. 



Figure 1-5. Sketch showing a typical scatter pattern. For most materials, the 
figure is rotationally symmetric about the slab. 
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Generally in radiative transfer work, one tries to expand the phase 
function in a series of Legendre polynomials (see the expansion of 
eq. (2-31)): 


N 

P {cos do) = ^ Cjj Py(cos^o) (u)o = 1) (1“34) 

j=0 

with cos 00 given by equation (1-33). It is quite obvious that the more 
forward scattering we have, the more terms in equation (1-34) may be 
required to accurately describe the phase function; i.e., N may have to 
be several hundred. 

In all of what follows, we shall not be quite so ambitious in our 
expansions. Practically all of the authorities from whose work most of 
the remaining text is drawn content themselves with at most two or 
three terms of equation (1-34). This makes the mathematics tenable 
and makes the physics of the process much more transparent in the 
resulting equations. Also, somewhat surprisingly, the numerical results 
are not too bad, and are of acceptable accuracy for many applications — 
for example, in climate modeling. 

Start with the one- term expansion 

P(cos0o) = l (1-35) 

which is obviously the simplest possible case, and which describes the 
very important case of isotropic scattering — i.e., scattering that is the 
same in all directions. The reader should not dismiss this simple case as 
being too elementary to be useful. Many radiative transfer processes are 
in fact very nearly isotropic and can be adequately studied by means 
of this analysis. Moreover, the comparatively simple solutions which 
follow from this assumption can be extrapolated to more complex cases, 
as the use of so-called “similarity” transformations frequently permits 
a transformation of variables from a more complex anisotropic case to 
an equivalent isotropic form. (See chap. 7.) 

Chandrasekhar presents some interesting results for the two-term 
expansion 

P(cos0q) = 1 + COS0O (1-36) 

The three-term expansion 

P(cos 0o) = 1 + cos 00 + 0 ) 2^2 (cos 0o) (1-37) 

is also of particular interest, as for the special case in which d)i = 0 and 
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0^2 = 1/2; this reduces to the well-known Rayleigh phase function, 

P(cos 0o) = ^ ( H- cos^ 0o) ( 

This phase function has equal forward and backward scatter peaks and 
is used to describe scatter phenomena by particles which are very small 
compared to the wavelength of the incident radiation. (See fig. 1-6.) 



Figure 1-6. Sketch of the Rayleigh phase function. The figure is rotationally 
symmetric about the long axis. 

Finally, the Henyey-Greenstein phase function is frequently used 
when a large forward-scattered peak is desired. This phase function is 
given by 


P(cos^o) = , 2 \ ^ — a \ 3/2 (1-39) 

(1 + ^/2 -2gcos6oy/^ 

where g is known as the asymmetry parameter, and controls the size of 
the forward peak. Equation (1-39) is particularly useful in theoretical 
studies involving asymmetric scattering because it is a generating 
function for Legendre polynomials and has the simple expansion 

oo 

P(cos^o) = 51 (^*^ + l)?”^n(cos0o) (1-40) 

n=0 


Positive g gives a forward-scattering peak, while negative g gives a 
larger backward- scattered component. In order to achieve a better 
approximation to a given phase function, two or more expressions of 
the same type as equation (1-39) or equation (1-40), with different 
values of g, could be combined. Note that g = 0 in equation (1-39) or 
equation (1-40) gives the isotropic phase function. 
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The ratio of the size of the forward peak to the backward peak can 
be found from equation (1-39) 


PjOp = 0 ) ^ fl+iY 

P{9o = n) \l-g) 


(1-41) 


from which table 1-2 can be extracted. Many aerosols have ratios of 
forward- to backward-scattering peaks of the order of several hundreds, 
and can thus be adequately represented for many purposes by the 
Henyey-Greenstein phase function with g of the order of 0.65 to 0.70. 


TABLE 1-2. RATIO OF FORWARD- TO BACKWARD- SCATTERING 
PEAKS FROM THE HENYEY-GREENSTEIN PHASE FUNCTION 


9 

P{9o = 0)/P(<?o = 7 t) 


1.000 


1.826 


3.375 


6.405 


12.704 


27.000 


64.000 


181.963 


729-000 


6859.000 


oo 


The azimuthal integration of equation (1-34) gives some particularly 
useful results. From the complete expansion given later in equa- 
tion (2-32), this results in 


1 

^(/^> ^ P{cos9o) d(t> (1-42) 

It can be seen from equation (2-32) that all terms except those for 
m = 0 integrate to zero over the range of 0 to 27T, and we are left with 

oo 

•P(m, m') ^ XI (1-43) 

y=o 
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and this in the two-term expansion gives 

/z') = 1 + Wi/i/i' (1-44) 

while the three-term Rayleigh expansion gives 

= 1 - 1) (1-45) 

a particularly simple and useful form. 

The (I?! in equation (1-44) is related to the asymmetry parameter, 

1 

(cos^o) ^ “ / P {cos 6q) cos Ood cos 6q (1-46) 

2 J~i 

for which, in the case of the Henyey-Greenstein phase function, 

(cos 9q) = g (1-47) 

For this case, 

= 3g (1-48) 

It is important to grasp the conceptual differences between scatter- 
ing and absorption. In the scattering process, the photon interacts with 
a particle of the medium in such a way that, macroscopically speaking, 
the direction of travel of the photon is altered, but (in all cases consid- 
ered in these notes) its energy remains constant. It can be imagined 
that the photon “bounces off” the particle in a particular direction, 
with no exchange of energy with the scatterer. Thus, neither the inter- 
nal nor the kinetic energy of the particle is changed, and consequently 
the “temperature” of the medium is unaffected by pure scattering. 

In the absorption process, on the other hand, the energy of the 
photon is completely transferred to the particle, and the photon ceases 
to exist in its original form. The kinetic energy of the particle is thereby 
raised — the “temperature” of the medium increases. Emission is the 
opposite of absorption. The medium particle ejects a photon and the 
particle loses energy — the “temperature” of the medium decreases. 

In general, a medium can absorb and emit radiation, and can scatter 
radiation, but only the absorbed or emitted portion of this energy, 
gained or lost from a given beam of radiance, can contribute to the 
energy change of the medium. In the present text, the term conservative 
scattering will refer to the process of pure scatter with no absorption 
or emission. 
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The Equation of Transfer 

We now derive the integro-differential equation which describes the 
total change in the spectral intensity, or radiance, as it traverses an 
infinitesimal distance through an optically active medium which can 
absorb, emit, and scatter electromagnetic energy in the wavelength 
interval du centered about i'. The equation will be derived first in 
a very general form, and then specialized to the various forms usually 
seen in the applications literature. 

Consider an absorbing, scattering, and emitting medium whose 
optical properties are characterized by a spectral volumetric absorption 
coefficient, Ki/(s), and a spectral volumetric scattering coefficient, 
where 5 is the distance along the absorbing path. A beam of 
monochromatic radiation of spectral intensity /(s, n, t) travels through 
the medium in the direction fl along the path ds, and is confined to 
the solid angle dfl centered about the direction fl. (See fig. 2-1.) We 
C8tn write the outgoing intensity as 

7^(s,n,0 4-D/^(5,n,0 (2-1) 

I^(s, t) ^ ^ I^(s + ds, t + dt) 

\J " 

dA ds 

Figure 2-1. The change in intensity of a monochromatic beam of radiation as it 

passes through an optically active medium of length ds. 

where the total differential term represents the difference between the 
intensity entering the left face of the slab and that leaving the right 
face. 

Let Wl/ denote the net gain or loss of radiation by the beam in this 
volume element per unit volume, time, etc. Then quantity 

Wu dA ds dU du dt 

PRECiJ>iNQ PAGE BLANK NOT FUi/!£D 


( 2 - 2 ) 
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represents the net gain of radiant energy by the volume element. But, 
by definition of the radiance, /j/, this is precisely equal to 


DIiy{sy n, t) dA dQ du dt 


(2-3) 


and hence 




= W^ 


(2-4) 


We have taken here an Eulerian approach to equation (2-4); i.e., we 
have assumed that we are stationary and are describing what goes on 
inside a fixed volume element dAds^hence, the use of the total or 
substantive derivation in equation (2-4). Equation (2-4) can be written 
in terms of the more common time and space derivatives by using the 
usual transformation 


Ds c Dt c dt ^ ^ 


where c = cfl is the velocity of light (the velocity of the photons). Thus 
equation (2-4) becomes 

I djAs^ n,t) ^ ^ ^ ^ 2 _ 5 ) 

The second term is simply the directional derivative of in the 
direction s, so we get 

C dt + 

The net energy gain, Wj/, can be broken down into four separate 
pieces: 

Wt/j energy emitted by the volume element 

Wi ^2 energy absorbed by the element 

W 1/3 energy scattered out of the volume element 

Wu 4 energy scattered into the volume element from all 
directions 


For now, let us simply denote the total energy emitted by the volume 
element into the direction fl by 




(2-7) 
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The contribution is given by equation (1-20), written per unit 
volume and solid angle 


= (2-8) 

Wiy^ is the loss of radiant energy scattered out of the incoming beam 
by the scatterers in the medium (see eq. (1-25)) 

W^^=-a^{s)Iu{s,VL,t) (2-9) 

and energy scattered into the beam, is given by equa- 

tion (1-30) 

Wu^ = ~au{s)j P{coseo)i^{s,n',t) dn' (2-10) 

O' 


Substitute equations (2-7) through (2-10) into equation (2-6) 


1 dIi,(s,Cl,t) 
c dt 


dliy(s, n,i) 
ds 


= - Ki,(s)ii,(s,n,t) - (Tu(s)iu(s,n,t} 


1 

47T 



P(cos^o)/i.(5,n',0 dQ' (2-11) 


and equation (2-11) is the radiative transfer equation (RTE) in its most 
general form for our purposes. 

For practically all atmospheric propagation problems, the first term 
on the left-hand side of equation (2-11) is many orders of magnitude 
smaller than the other terms, and can safely be dropped from further 
discussion. 

The term which represents energy added to the emerging 

beam by emission, and the integral scattering term, which represents 
energy added to the emerging beam through scattering, are usually 
combined to give what is usually referred to as the source term, 
which represents the total energy added to the beam by emission and 
in-scattering: 


3 u = ju(s,t) + ^au{s) j P{cos6o)Ii,{s,{l',t) dQ' (2-12) 
n' 
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and hence, we can write equation (2-11) in the form 


dii/{s, n) 

ds 




(2-13) 


in which the time-dependent term has been dropped. Note that we have 
switched to a total derivative notation in equation (2-13). This is not 
strictly correct, as /j/ is a function of more than one variable. However, 
this is a convention that has been adopted in the RTE literature, and 
hence will be adopted here. 

Finally, we divide through Ki/{s) and write equation (2-13) 


as 


1 dI^{s,Q) 
Ku{s)+cru{s) ds 


+ /^(s,n) = Jt/(5,n) 


(2-14) 


where 

^ Ku{s)\a^{s) 

is referred to as the source function. 

Equation (2-14) is very general. We now make a very important 
assumption, namely, that the volume element is in local thermodynamic 
equilibrium (LTE) with the surrounding medium. This LTE assumption 
is valid in most atmospheric problems, at least below 30 to 50 km. Then, 
Kirchhoff’s law allows us to define in terms of the Planck function, 
Bu{T), 

jt = K^{s)Bu{T) (2-16) 


where T is the absolute temperature of the medium in the volume 
element dAds. The source term then becomes 


ju = K^{s)B^{T) + ^a^{s)f P{cosOo)I^{s,n') dn' (2-17) 

n' 

Define the spectral volumetric extinction coefficient^ 

P^{a) = Kt,{s) + a^{s) (2-18) 


and the ratio 


= 


<^i/(a) 


(2-19) 


Ui/ is called the single-scattering albedo^ or the particle albedo, and 
expresses the fraction of the attenuated beam which is lost to scattering 
alone. In terms of cDi/, 
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Figure 2-2. Geometry of plane-parallel atmosphere. The direction z is mea- 
sured upward from the planet surface. Positive p denotes upward-traveling 
radiation. 

and hence, we can write equation (2-14) as 

+ = ( 2 - 21 ) 

where the source function is 

Ms,n) = [\-uu)B[T) + p[cos 9 )iu{s,n') <m' ( 2 - 22 ) 

n' 


The RTE in Plane-Parallel Atmospheres 

We now further confine our attention to the passage of radiation 
through a plane-parallel medium. This is a medium which is stratified 
in planes perpendicular to a given direction such that the optical 
properties of the medium are functions of and u only. Since the 
thickness of a planetary atmosphere is generally small compared with its 
radius, this assumption is almost universally made in applications of the 
RTE to atmospheric radiation studies. Now, d{ )/ds — cos 6 d{ )/dz = 
fxd{ )/dz (see fig. 2-2), so that we can rewrite equations (2-21) and 
(2-22) in terms of /z, and 0: 

^ 4 ,) = (2-23) 
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Here, in 

dn' = sin e' de' d<t>' ^ -dp! (2-24) 

so that 


0 ) = (1 - Cji,)Bi/[T{z)] 



P(cos$Q)Iu{Zy d<t>' 


(2-25) 


It is also convenient at this time to introduce the concept of optical 
depths Tiy, defined to be 


Ti, = 



dz^ 


' dxiy = —/3i/{z) dz (2-26) 


Note that the optical depth is defined to be zero at the top of the 
atmosphere, and increases as one descends through the atmosphere, 
in the opposite direction from that in which z is defined. This 
convention is a carry-over from the astrophysical literature, where, 
in studying the radiative properties of stellar atmospheres, distance 
and optical parameters are measured positive from the surface of 
the star inwards. Since much of radiative transfer theory has been 
developed and published in connection with studies of stellar interiors, 
this convention has, for the most part, been adhered to in applications 
of radiative transfer theory to planetary atmospheres. 

If the height variable z is replaced with the optical depth Tt/, in 
equation (2-23) (see fig. 2-3), then 




diu {Tu,fJ;4>) 
dTu 


"H Ii/ (^i/, /i, <t>) — Ju 0) 


or 




dll, 

dri. 


Ii/ {ti,, /i, (f>) — Ji, (tv, /i, <f>) 


(2-27) 


with 


Ju <^) — (1 ~ ^^1^) [T { ti /)] 


+ 



P(cos0o)7i/ (riy,^',0') dp d(j>^ 


(2-28) 
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Figure 2-3. Sketch showing the relationship between the vertical coordinate z 
and the optical depth, Tu- 


Further Specialized Forms 

Two further specializations of equations (2-27) and (2-28) are fre- 
quently encountered. 

First is the case where the emission term, Bi/ [T (r^^)], is small and 
can be neglected. In this case, only scattering and absorption are 
included in the transfer process; this is the situation usually encountered 
in studying radiation emitted directly by the Sun. This radiation is 
absorbed and scattered by the Earth’s atmosphere, but the atmosphere 
itself is cold compared to the Sun, and thus, at solar wavelengths its 
radiation is small compared with that emitted by the Sun. Thus, when 
the emission term is small, equations (2-27) and (2-28) are usually 
written as the single equation 


d/t/ 4 >) 

^ dTu 

r\ 

= 7i/ ^ J j P(co^eQ)Iu dy! d4>^ (2-29) 

Equation (2-29) and sundry of its equivalent forms will be the 
starting equation in much of what follows in these notes. 

The second specialized case for equations (2-27) and (2-28) occurs in 
the IR spectral region, where it is the scattering which can be neglected 
(except in clouds). In the case of measuring IR radiation from the 
Earth’s atmosphere, the emission term is the only source of radiation, 
and hence must be included. For this case, Cju = 0, and equations (2-27) 
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and (2-28) become 




dll/ ('7’t/j 0) 

dri, 




(2-30) 


Expansion of the RTE into Legendre Polynomials 

Equation (2-29) is still extremely difficult to solve. Part of this 
difficulty is due to the azimuthal dependence of It/ through the phase 
function. By expanding the phase function in a Legendre polynomial 
series, the azimuthally dependent terms in the function can be uncou- 
pled. The form of the expansion will then suggest that the radiance 
should be expanded as a Fourier cosine series. The result of substi- 
tuting these expansions in equation (2-29) is a set of uncoupled linear 
integro-differential equations for the various orders of expansion. From 
this, we will show that only the azimuthally independent equation con- 
tributes to the flux calculations. Since this is the parameter of great- 
est interest in most atmospheric applications, we can then confine our 
future attention to the solution of only this azimuthally independent 
equation. 

First, we expand the phase function, equation (1-29), in a Legendre 
polynomial series of order N: 

N 

P{cos9q) = ^ Qj Pj (cosBo) (u)o = 1) (2-31) 

where cos0q is given by equation (1-33). Then by the addition theorem 
for Legendre polynomials we can write equation (2-31) in terms of 
/i, pf <j>, and (f)^: 


N N 

P(/z,/z',(;6,0') = E E"f cos[m(^' - 4>)] (2-32) 

m=0i=m 


where 


wf = (2 - <5^) 


{£ — m)! 
{£ + m)! 


/ 0 < m < N 

\£ = m,m + I, . . . , N 


) 
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Substitute equation (2-32) into equation (2-29) 


N N 


m=0 

n l 

cos[m(0' - <^)] d/i' dd>' (2-33) 


Note that the phase function has separated into the product of a 
function of fi and only, times a function of {(f>^ — (p) in each term. 
This suggests that if we expand 0) in a Fourier cosine series 

in we ought to be able to separate the azimuthally dependent terms 
from the azimuthally independent terms by equating like coefficients 
of cosm(0' -- 0). The direct sunlight, which is usually taken to be the 
source of the radiation in the atmosphere, is assumed to be directionally 
defined by the angles 0 o) (see fig. 2-4). Since most of the radiation 
will be along this direction, let us expand about this unit vector 


N 

cos[m (</)0 - (t>)] (2-34) 

m=0 


where the coefficient IJP is a function only of Ti, and /x but not of <j>. 



Figure 2-4. Sketch showing the scattering of an incoming collimated beam of 
solar radiation. 
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Substitute equation (2-34) into equation (2-33) 

> , cos[m(0o - ^)] 

m—O 

N ^ N N 

m=0 m=0 i=m 

N 

X cos[p{(f}Q - (j/)] cos[m(0' - (f>)] dp! d<t>^ (2-35) 

p=0 



N .1 

E/ 

P=0''“^ 


Now 


Examine the integral term of equation (2-35) 

r27T rl ^ 

/ / E ^J(^*^’/^')cos[p(0o - 0')]cos[m(0' - 0)] dp* d<f>* 

Jq J-\ 

r2TT 

dp' I cos[p( 0 o - 0 O)cos[m{<^' - <^)] d(p' 

Jo 

(2-36) 

r2n 

/ cos[p(0o ~ 00] cos[m(0' — 0)] d(f)^ 

Jo 

= 27t (p = m = 0) 

= ncosm{(f>Q — (f>) (p — m 7 ^ 0) 

= 0 {p^rn) 

Thus, we are able to write the right-hand side of equation (2-36) as 

(1 -h ^0^) 7TCosm(0o “ j (2-37) 


If we now substitute equation (2-37) into equation (2-35) and equate 
coefficients of cos[m(0o ~ 0)1 on both sides of the equation, we can write 
for the IJP component 


dljf[ru,p) dp 
dru 


N 


= C(r.,M) - X (1 + ^"*) E 


/: 




(2-38) 
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The definition of spectral fiux, equation (1-14), can be written as 
r27T rn 

= / / Iu{Tu^&y<f>)<^os6smO dO d<f) 

Jo Jo 


or 


r27T r — I r '^TT r I 

Fi,{Tiy) = - / pJ ,(!>) d(f) = I / ,<i>)dy! d<^> 

Jo Ji Jo j-l 

(2-39) 

Substitute equation (2-34) into equation (2-39) 


r27T r 1 


fJ' yZ /^(T-i/,//)cos[m(0o - <t>)] dpi d(j) 

-1 r^O 

ri r2TX 

= y ^ ptlJ^irjy.pi) dfjL cos[m ((^0 ~ 0)] # 


But this vanishes unless m — 0, in which case we get 

= 2n J ^ pt) dpt (2-40) 

This demonstrates that the flux depends only on the m = 0 term— that 
is, the azimuthally independent term of equation (2-38). So, we will 
now restrict all future developments to equation (2-38) with m = 0 and 
drop all the superscripts: 


^ e^o 

or in a somewhat prettier form 

J ^ Iu{Tu,p')P{p-,p!) dp' (2-41) 

where 

N 

'P(m> /) = XI (2-42) 

e=Q 
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Equations (2-41) and (2-42) are the forms most frequently seen in 
the literature. Remember the restrictions, however: 

1. No thermal emission, 

2. plane-parallel atmosphere, 

3. phase function expandable in Legendre polynomial series, and 

4. azimuthal symmetry. 


RTE for Diffuse Component Only 

We now derive one other form frequently seen in the literature. In 
the preceding development, the term lu was considered to be the total 
spectral intensity. In problems of atmospheric physics, the assumption 
is usually made that the Sun’s rays consist of a parallel, or collimated, 
beam of radiation hitting the top of the atmosphere at some direction 
specified by the angles 6o and 0o* Some of this radiation is multiply 
scattered and appears at various values of Ti, in the form of diffuse 
radiation; i.e., radiation which has been multiply scattered and is now 
traveling in all directions. Another portion of the incoming solar beam 
is absorbed by the intervening atmosphere between the entry point and 
the current value of The remainder appears at as attenuated 
solar radiation. This component is referred to as the direct component, 
traveling in the same direction as the incoming beam. In analysis, 
it is frequently convenient to separate these two components in equa- 
tion (2-29) so that the resulting equations describe the behavior of 
the diffuse component only. This also simplifies in many ways the 
application of the boundary conditions. 

So, let us write 

/^ = /^ + jS (diffuse + solar) (2-43) 


As indicated above, the direct beam consists of photons which were 
originally in the incoming solar beam. These represent what is left over 
after all the scattering and absorption has taken place. It does not 
include photons which have been scattered out of the incoming beam 
and then scattered back into the original direction — these are part of 
the diffuse component. It also does not include photons emitted by the 
layers of the atmosphere above it— these are also part of the diffuse 
beam when, rarely, the emission terms are included in the RTE. Thus, 
the direct, or solar, beam satisfies its own differential equation, of the 
form 


4 >) 

^ d,Tu 




(2-44a) 
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with the upper boundary condition prescribing the incident radiance to 
be a beam collimated in the direction (0Oi </>o) 

/^(O, -po, 4 >q) - t^Fq8{p - Mo) S{(j) - <j>o) 

Solving and applying the boundary condition yields the intensity for 
the direct beam in terms of the incoming solar flux, n F q 


-p, <i>) = TTFoe-^‘'/f^6{p - po) 6{<f> - <f>o) (2-44b) 


where the 6 are Dirac delta functions. The factor tt is frequently 
introduced into the solar flux because of the way Chandrasekhar defines 
the flux term. He defines the flux as 





pIu{tu,p) dp d<j) 


rather than our definition in equation (2-39). The reason for this is that 
the factor tt then cancels out of both sides of many of the flux equation 
forms, thus eliminating the necessity of carrying the 7r-factor through 
a lot of theoretical development. 

Now, put equation (2-43) into equation (2-29) 


diu (tu,u,<I)) dl^ ( tu,p,(1>) jD, , tS, .. 

CLTif dTif 

f2TT rl - 

+ dfi' d<j}* (2-45) 

From the differential equation (2-44a), the second term on the right- 
hand side and the second term on the left-hand side of equation (2-45) 
are equal. Substitute equation (2-44b) into the part of the integral 
term of equation (2-45) 

/'27T rl 

/ P{p,<l>-,P,<t>')^Foe~"-'^^°6iP -PoW -Mdp' d<j>' 

Jo J-i 

= (2-46) 
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Use equation (2-46) to write equation (2-45) in the following form, 
dropping the superscript D, where I„ is now understood to be the 
diffuse component only 

^ ~ ~ J J ,<t>) dp' d<j>' 

-^Foe-^''!'^<^P{p,<l>--po,4>o) (2-47) 

Note that for the special case of isotropic scattering, P{p, (j>; p' , 0') = 1, 
dlu{ru,p,4>) _j j ,, f f , , ' i'\ 1 ’ j'> 

^ ^ y'y4>) J J ) dfi d(j) 

- '^Foe~^''/'‘°P[p,4>-, -po,<l>o) (2-48) 

and if in addition we assume azimuthal symmetry 

^ u)i/ T f\ 1 f 

^ 2 J ^ /^ ) 

- ^Foe-^-'/f^Pip, -po) (2-49) 

The azimuthally symmetric form of equation (2-47) is 

dli/{ji/^p'^ Jf f !\ T f f\ 1 f 

^ d^ - Y J ^PiF,t^)Iu{Tu,p) dp 

-'^FQe-^^l^'^P{p,-^Mf) (2-50) 

Substituting equation (2-42) in equation (2-50) gives 

^ N I 

^djA^±2Fl f Pf{/j,')I^(T:y,p') dp' 

" t=Q •'-1 

N 

- (2-51) 

This equation is frequently used as the starting point for the devel- 

opment of the diffuse components of the two-stream and Eddington 
solutions to be derived in chapters 5 and 6. 
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Notice the obvious differences between the forms for the azimuthally 
independent equations, (2-41) and (2-50), The latter equation contains 
the solar flux exponential term while the former does not. The presence 
of this exponential term is a giveaway that the radiation term contains 
the diffuse component only, whereas the form of equation (2-41) includes 
both the direct and the diffuse components. This difference, while 
obvious now, is not always pointed out by authors in the open literature 
and, hence, a careless application of their equations may result in rather 
strange-looking results, especially when applying boundary conditions. 


37 




Chapter 3 


Formal Solutions to the Intensity and Flux Equations 

Return now to the form of the RTE in a plane-parallel atmosphere, 
equation (2-27) 


(3-1) 

uTj/ 

and derive the formal solution to this equation. It should be remarked 
that this is not really a solution to equation (3-1) in the normal sense 
of being used to derive very complex numerical results. It can, in fact, 
be used in some very simple cases, as will be demonstrated later, but 
in general the coupling between the intensity and the source function 
precludes the formal extraction process of yielding 1^ as an explicit 
function of and fi. The formal solution’s utility is that it forms the 
starting point, either in the spectral intensity form or the flux form, 
both for theoretical analyses and for some elementary methods. 

It is convenient to break the solution into two parts, one for the 
upward component 

(o</i<i) 


and one for the downward component 

(-l</i<0) 

with the solution subject to the boundary conditions 


7^(0, -fi, </)) = ll{0, -fj) 
(^1/ 1 Ml *A) 


PRECEWNQ PAGE BLANK NOT 


(3-2) 

(3-3) 
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Note now that 


A. 

dru 




so if we divide equation (3-1) by p and multiply by e equa- 

tion (3-1) becomes 


Let us integrate this equation between the general limits t\ and ^2 

Jtl ^ 

Now, we want to find the upward component of radiation at Ti^. 
This radiation comes in part from the surface at and also from 

all of the infinitesimal layers of the atmosphere between and r*, 
all properly attenuated by the intervening layers of atmosphere, (See 
fig. 3-1.) So, in equation (3-4) we let and t 2 = r* and solve for 

^i/('?Vj 0) 


Jtu M 


(3-5) 



] 


//////////////////////////A 

277777777777777777777777777 




V 

Figure 3-1. Upward radiation at Tu due to radiation from the surface at r* and 
from intermediate layers of atmosphere at r^. 

Similarly, the downward component of intensity at r^/ is equal to the 
downward intensity impinging on the upper boundary at = 0 plus 
all the source terms between the top layer and Tu>, also attenuated by 
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the atmosphere. (See fig. 3-2.) Thus, in equation (3-4) we let <i = 0 
and <2 = and solve for p, (f>) 

/i;(T-iy,Ai,(A) = lU^,p,<P)e^''^^ 

- r 4>)^ (3-6) 

Jo 

in which —1 < ju < 0. 






/innimumnnnnnnnu/rn t ^* 


Figure 3-2. Downward radiation at Tu due to radiation impinging on the top 
surface = 0 and the intermediate layers of atmosphere at ri- 


Equations (3-5) and (3-6) demonstrate the usual exponential nature 
of the attenuation of monochromatic radiation with increasing optical 
depth. This requires, of course, a negative argument in the exponential, 
while equation (3-6) appears to produce a positive argument. However, 
p is negative, thus giving the proper sign. So in order to make 
the equation look right, most atmospheric physicists at this point 
replace the p with —p and incorporate the minus sign explicitly in the 
exponential term of equation (3-6). This convention tends to complicate 
the interpretation of the ensuing equations to some degree, but as it is 
a relatively minor nuisance, and is consistently done in the literature, it 
will be followed here also, with appropriate caveats as the need arises. 
Then equation (3-6) becomes 

-h r J.(r^, -p. <j>)^ (3-7) 

Jo M 

Equation (3-7) is the desired equation for the downward component of 
spectral intensity. 
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To get the flux components, we proceed from equation (2-39) 

r27T r\ 

Fu{tu) = y dn d4> (3-8) 

which we break into components as 



At this point we must be wary in the literature. If the convention 
— adhered to, then equation (3-9) can be continued 

directly as 


FAtu) = Fj(r^) Fi(r^) 

(3-10) 

with 


r2n rl 

dp d(f> 

JQ Jo 

(3-11) 



2 0 

d(j> 

jQ J-1 

(3-12) 

But if the convention 0 < /i < -fl for // is followed, and we replace 
p with —p for the downward component, then we proceed from equa- 
tion (3-9) as 

F,y{Tu) = [ f dn d(j)+ f f 

Jo Jo Jo J -1 

Py (f>) dp d(f> 

= / [ dfi d(t>- f f dfi d(j) 

Jo Jo Jo Jo 

= [ f dn d<l>- f f 

Jo Jo Jo Jo 

—py (t>) dp d<}> 

= FI{t^) - F^{t„) 

(3-13) 

where ^ 


^J(7v)=/ / dfi d4> 

Jo Jo 

(3-14) 

and ^ ^ 


FH'^i') = / dfi d(j> 

Jo Jo 

(3-15) 
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Note the difference in signs and the difference in the integration limits 
between equations (3-10) and (3-13). is evaluated in the same 

way in both cases, but is handled somewhat differently. So, we 

will develop fI{ti,) first, and then develop both expressions for F^{tu)- 

fI{ti^) is given by either equation (3-11) or equation (3-14), using 
from equation (3-5) 

^27T ^ 

dp d(j> 

Jo Jo 

f [ ^ [ dp d4> (3-16) 

Jo Jo Jtu ^ 

This is about as far as we can go analytically with equation (3-16) 
without having any knowledge about either the directional distribution 
of the radiance or the source function. We can proceed with this 
development if we make the assumption that the phase function, and 
hence the source function, are isotropic. This is a rather limiting 
restriction, and applies only to the case where the source function 
can be replaced by the Planck function, as in equation (2-30). This 
then becomes a problem of emission and not scattering. The resulting 
equations are not applicable to general scattering problems, but are 
applicable to studies of infrared radiation. These forms are frequently 
seen in the literature, and it was thought not unreasonable to present 
them here, even though the main thrust of these notes is with the 
scattering problem. 

In any event, we make the above assumption and write equa- 
tion (3-16) as 

= 2'kII{t*) f dfi 

Jo 

+ 2it f drl f dfi (3-17) 

Jti/ Jo 


The /x-integrals are exponential integrals of various orders, where 





e 


—xt 




dt 


(3-18) 
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Evaluation of the first integral: 
Let 

1 


di 

/i = 7 = 


then 


Jo Joo ^ ^ 


fOO 

= ^3 = Esir^ - T„) 


Evaluation of the second integral: 
Let 


f dfi = - f' 

Jo Joo ^ 


-L 


oo 


d^ = E2{rl~T^) 


so that we can write: 


FI{tu) = 2nllT*E3{T* - r„) + 2-ir E 2 {tI - t^) Jt/(r') dr^ (3-19) 

J Ty 

This equation could stand as it is, but it is more convenient to put it 
into another form — from equation (3-18) 

dEn{x) f°^ e 


dx 


too roo 


Thus 


dEzirl - T^) 


drl 


'dEsirl - TuY 

■(/(r/, - r^)' 

. d{Tl, - Tu) . 

dr!. 


-E2{t'u - Tu) 


so that we can write equation (3-19) as 

Fj(r^) = 27r/J(r;)£3(T* - t ^) 


2-k 


/ Mrl) 

J Ty 


UTj, 
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Integrate by parts and write equation (3-20) in the form 


fI{tu) =27rE3(r; - r^) [/^(r;) - JuK)] + 27rJ^M 

+ 27t E^{tI - (3-21) 

The reason for bringing in the £*3 in the second integral rather than 
leaving the form as £2 Is that, as we will see later in equation (4-26), 
2E^{tu) is an angular integrated monochromatic transmission function, 
Tr(T|y), which, when integrated over frequency, can in some cases easily 
be evaluated from band transmission models. Guided by this concept 
then, we write equation (3-21) in final form as 

= 7rTr(r* - Tu) + 27rJi/(Tt,) 

+ 7T fr{rl - (3-22) 

Note the physical difference between the two terms in the bracket: 

nll{r*) == flux from the surface at Tjy = r* 

= source flux from the atmosphere immediately 
above the surface at Tj/ = 


Now we evaluate the downward flux components. This is done 
exactly as above for £j(v), except that there are some mildly tricky 
steps involved in the manipulation of the £-integrals that can easily 
give wrong signs if one is not very careful. 

First, we will evaluate equation (3-12). Then we must use the 
defined by equation (3-6) 


n o r rTu ^ 

-1 L Jo 




dfj. d4> 

(3-23) 

and if we make the isotropic assumption on and Ju this reduces to 

/ O rTu rO ^ 

dn-2w / MtI) drl / dfi 

(3-24) 
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Evaluate the first integral: 



Let 




dti = -^ 


Evaluate the //-dependent part of the second integral: 


Again convert E 2 to an E:^ derivative: 

dE^[rv - rl) _ dE^(ri/ - tI) d{ry - tD 


drl 


d[jv - tI) drl 


= £^2(riy - tI) 


Hence, equation (3-24) becomes 

= -2-rll{Q)Ez{Tu) -2n f drl 

Jo 

Integrate by parts and write equation (3-25) as 


Fj(r^) = - ~ + 2T^Ju{(S)E^[r^) 

+ 27T r Esir,-Ti)^^^0^drl 
Jo drl 


or, in terms of the transmission function 


= - ttTVCv) [4(0) - J^(0)] - 2i^Mtu) 


+ 


fTv ^ 

7T / Tr{T„ 

Jo 


) \ dJi/[ji/) j I 


<) 


dr!, 


(3-25) 


(3-26) 


(3-27) 
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The net flux follows from equation (3-10), using equations (3-22) 
and (3-27) 

=nfr{T* - v) [llK) - J^(r;)] 

- TrfM [/^(O) - MO)] + 7V(r^ - 

dJi>(ju) j / 


+ 7T y lr{Tjy — Ti/) dTf^ 


Now, 7 V(t^ - M = Tr{rl - r^), since these are transmission functions 
between the optical depths and r' , and are assumed to be the same 
numerical value when taken in either direction. Thus, we get 

FMi') ='^Tr{r* - M [ 4 (r*) - Ji/(r*)j 
- TrfriM [/i(0) - J^(0)] 


+ 7 rJ^ r,(v O M 


(3-28) 


This equation for the net flux forms the starting point for many studies 
of the temperature structure of the Earth’s atmosphere, and is used to 
describe the infrared cooling part of this structure. See, for example, 
Rodgers and Walshaw (1966). 

Now, we develop from equation (3-15), where in this case we 

must use equation (3-7) to define —fi, 4>) 


n l 

I 

= 27r/i(0) f 

Jo 


/J dfi 


Jo ^ 


d/i d(f> 


fie 


-Tu/fl 


d^ + 2n r MrDdri rf// (3-29) 

Jo Jo 


Evaluate the first integral: 


= -1/^ 


/„ 


1 


fie dfi — E3{Tiy) 
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Evaluate the /i-dependent part of the second integral: 

f dn = E2{T^-Tl) 

J 0 

so that equation (3-29) becomes 

Fi (r,) = 27tII{0)Ez{t,) + 27T r drl (3-30) 

Jo drl 

Integrate by parts and equation (3-30) becomes 

Fj(r^) =27r£'3(r^) [/i(0) - J^(0)] -h 27rJ^(r^) 

- 27T E^{t, - rl)^^0^drl (3-31) 

or 

= T^fr{Ti,) [/J;(0) - Ji/(0)j + 2-nJ^{Ti,) 

jj Tr{T^ - (3-32) 

It can be seen that equation (3-32) is exactly the negative of equa- 
tion (3-27) which is extremely fortunate or we would have a serious 
problem in computing the net flux — and hence, the net flux, given now 
by equation (3-13), also produces equation (3-28). 

The message here is that, when reading the literature, one’s atten- 
tion must be drawn to the way the author defines the downward flux; 
i.e., whether with -1 < < 0 or with 0 < /x < -|-1, and, hence, whether 

the net flux is defined by equation (3-10) or equation (3-13). 
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Reflection and Transmission Coefficients, Surface Effects, 
and Albedo 

In many applications of radiative transfer theory, we are not partic- 
ularly interested in what goes on in the interior of the atmosphere, or 
inside a finite thickness of the atmosphere. For instance, we may be in- 
terested only in what comes out of the top of the atmosphere at = 0, 
or what comes out of the bottom at t^/ = In order to simplify the 
extraction of these data from the radiation field, a theoretical approach 
known as the principle of invariance was developed by Ambartsumyan 
(1958), and further developed and clarified by Chandrasekhar (1960). 

We mention this principle here in order to provide a springboard 
for introducing the concepts of reflection and transmission functions, 
to which this chapter is devoted. The principle of invariance will be 
discussed in detail in chapter 8. For now, we merely state that if 
the reflection and transmission properties of two thin slabs of optical 
material are known, then the principle allows us to determine the overall 
reflection and transmission properties of a composite slab made by 
placing the two thin slabs face to face. The solution to the radiative 
transfer equation for a thin slab is relatively simple (see chap. 5). 
Thus, when working with atmospheric problems we could divide the 
atmosphere into a number of thin layers, use the thin-layer solutions 
for the RTF to determine the transmission and reflection properties 
of the thin layers, and then use the principle of invariance to build 
up the atmosphere layer by layer, and thus compute the reflection 
and transmission properties for the finite-thickness atmosphere without 
having to solve the complete form of the RTF. This is a particularly 
useful concept in deriving numerical results for both homogeneous and 
nonhomogeneous atmospheres. 

We proceed now with the introduction of the transmission and 
reflection functions. 

Chandrasekhar defines the scattering function, 

S{t* : n,4>,fio,<l>o) 
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and the transmission function 

f(r* : 


by the following equations 


1 

/ / 5(r* : 0, /i', 0 ')/jNq{/, 0') d/ 

1 Z*^ - 

^TRANS(^*> J J T(t* : dfx^ d4>^ 


(4-1) 

(4-2) 


where /ref the reflected diffuse radiation, /trans f^he transmitted 
diffuse radiation, and /inc the incident radiation. Both S and T are 
explicit functions of the total optical depth, r*. The factor l/fi was 
introduced to secure the symmetry of S and T in the pairs of variables 
iH,4>) and i/iio,(/>o); i.e., 


S(t* : fi,</>,^o,^o) = S(t* : fiQ,<po,M,4>) 

T(t* : n, (p, (M), 4>o) - T{t* : fiQ, (pQ, fx, <p) 

If the incident radiation is considered to be solar radiation, entering 
the atmosphere in a parallel, or collimated, beam, then we can write 
(see eqs. (2-44)) 


^INc(/^^ 4>') = TrFo6{(x' - fio)S{(p' - <po) (4-3) 

where ttFq is the solar flux, and the 6 are Dirac delta functions. 
Substitution of equation (4-3) into equations (4-1) and (4-2) gives, for 
a collimated incident beam, 


4>) = : fi,(p,fJa,^) (4-4) 

Fq~ 

^TRANs(^*i “M. <l>) = ^0’ ^) H'^) 

Note that T defines the diffuse component of the transmission only. 
The reduced direct component, is not included in T. 

There is another set of definitions of reflection and transmission 
coefficients which appears frequently in the literature (e.g., Liou, 1980), 


50 


Chapter 4 


which is defined by the relations 


1 /’27T r 1 

7T Jo Jo 


(4-6) 

1 /’27T Z* 1 

^TRANS(^*> -Ml <^) = “ / / ^{m, <!> '■ M^ 0O^INC(-M^ <A0m^ 

(4-7) 

Note that these differ from S and T in that they drop the \/n in front of 
the integrals, but include a /i under the integral. Thus, the integrals in 
equations (4-6) and (4-7) are closely related to flux integrals. This is the 
reason for the I/tt factor, changing the integral from flux to radiance, 
as required by the left-hand side. 

If we again use equation (4-3) for the incoming flux, we get 


Jref( 0, Ml 4>) = <p : fiQ, (f>o) (4-8) 

-fTRANs(''"*i — Ml <A) = MO-^o'^(Mi 4> '■ MOi 4>o) (4-9) 

Comparison of equations (4-4) and (4-5) with equations (4-8) and (4-9) 
gives the correspondences between the two sets of coefficients 

(4-10) 

4/ioM 


T = 


T 

4/xqM 


(4-11) 


Frankly, it is not at all obvious which set, if either, is better to use — 
R and T, or S and f. There seem to be some small practical advantages 
in using R and T, since they are defined in terms of fluxes rather than 
radiances, and since the albedos are generally also defined in terms of 
flux quantities. The set S and f seems to be used more frequently 
in high-powered theoretical developments than does the set R and T, 
but this may be due more to the impetus given these parameters by 
Chandrasekhar’s studies and influence than to any inherent advantage 
of their own. 

From the definitions of equations (4-6) and (4-7) we can immedi- 
ately write expressions for the diffusely reflected flux and the diffusely 
transmitted flux 


^REf(0 


n l 

J 


R{fi,(p ■ M^0O^INc(-M^<^OM^ dfj,' d<t>' (4-12) 
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r2n rl 

^TRAmiT*, -M, 4>)= / T{n, (j> : fi', <?!»O^INc(-/. dfj! 

JO Jo 

(4-13) 

We define the planetary or local albedo as the ratio of the total 
outgoing flux at the top of the atmosphere to the total flux entering the 
atmosphere at the angles 0q and <i>o. The total incoming flux is 


^INC = 




dp^ d(j/ 


and for the collimated beam of equation (4-3) 


^INC = ttFoM (4-14) 

The total outgoing flux is found by integrating equation (4-8) over all 
angles p and 0: 


^REF ^ <!> * llo^ <t>o) dp d<j> (4-15) 

and hence, the planetary albedo, r(//o), is given by 

„C,, 'I _ '^REf(0,/X,<?^) ^ nf i /’I , ,, 1 --, 

^l/^o) ^ — ~\ / R{fi,<j> : dfi d(^ (4-17) 

^INC Jo Jo 

Similarly, the diffuse transmission function, t(/io)> is written 


= 


7^trans(^*i 

^INC 
1 /■! 

- / T{fi, <f > : Ho, <^)/x dH d(f> 

7T Jo Jo 


(4-18) 

(4-19) 


Again, for emphasis, equation (4-19) describes the diffuse transmission 
function only. The direct transmission function is given by e~'^* 
However, the more fundamental definition in equation (4-18) may 
include both transmission components. 

For the special case of azimuthal symmetry, equations (4-17) and 
(4-19) reduce to 

r{no)^2f R{n,Ho)tJ'dn (4-20) 

and 
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t{no) = 2 [ T(/z,/io)M dfi (4-21) 

Jo 

respectively. 

The spherical albedo is defined for a planetary atmosphere as the 
ratio of the total flux reflected at all angles by the planet to the total 
flux incident on the planet. If we let the radius of the planet be a, the 
total flux incident on it is 

(7rFo)(7ra2) (4-22) 

We want to find the total flux reflected by the planet. Let dA be the 
area of the elemental ring on the surface of the planet as shown in 
figure 4-1. 

dA = 27ra^ sin 6q dOo 



Figure 4-1. Sketch of the geometry involved in computing the spherical albedo. 
The elemental area normal to the incoming rays is 
dA cos Oq = 27ra^ sin 6q cos 6q dO^ 
or, with ^0 = cos Oq 

Pq dA = ~2a?po dpo 

Thus, the total flux reflected by the element dA is 
(-27ra^^O duoj k^b»‘(Ato)] 
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and integrating over po from 1 to 0 

^REF = (27ra^)7rFo / Mor(Mo) dfiQ (4-23) 

JO 

The spherical albedo, f, then becomes 

- ^ [2na?)TrFo fM)T{fio) d/iQ 
TKOL^ • TtFq 

f = 2 / por{po) dpo (4-24) 

Jo 

In a similar manner we can define a spherical transmission function 

t = 2 f pot{po) dpo (4-25) 

Jo 

The spherical transmission function for the direct component is 

to = 2 f <Jf^o = 2i?3(^*) (4-26) 

(See eq. (3-18).) We see that is related to the direct transmission 
and that this is the reason for changing from E 2 to £*3 in the develop- 
ment of the flux equations in the last chapter. 

Inclusion of Surface Effects 

The reflection functions and albedos we have derived so far are 
for the atmosphere alone. If the atmosphere is bounded below by a 
reflecting surface, as it obviously is, then the reflection function and the 
albedo of the total system must be modified somehow by the presence 
of the surface. We now consider this problem, and use the approach 
of Tanr^ (1982). This approach permits us to work the atmospheric 
problem alone, without considering the surface effects, and then add 
the surface effects separately. Thus, the optically thin atmosphere and 
the single-scattering solutions introduced in the next chapter acquire 
considerable importance. 

Consider an atmosphere of optical depth r* bounded below by 

a Lambertian surface; i.e., a surface which reflects equally in all 

directions. Assume that each point of the surface is Lambertian with 
a reflectance and that each part of the surface receives the same 

downward flux. The solar flux at the top of the atmosphere is, as 
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usual, denoted by ttFq, and it enters the atmosphere at the angles no 
and The total flux received by each surface element is the sum of 
three separate components (as shown in fig. (4-2): 

1. A direct flux component — the incoming solar flux attenuated along 
the slant path, 

7rnoFoe-^*/>^o (4-27) 

2. A diffuse transmission component, arriving at the surface after 
multiple scattering (see eq. (4-19)) 

7T/xo^O^{-Mo) (4-28) 

3. A diffuse component arriving at the surface after multiple scatter- 
ings and reflections between the atmosphere and the surface. 



Figure 4-2. Sketch illustrating the three ways a specific photon can interact with 
the surface. 


Write the total transmission through the atmosphere as 

Tri-no) = + (4-29) 


Then, the total flux which reaches the surface before any surface 
reflection occurs is 

TTfioFoTri-no) ( 4 - 30 ) 

and the total flux reaching the surface after multiple reflections and 
scatterings between the surface and the atmosphere is 


7TfioFoTri-no) 


Psf + + plf^ ■!-•• •] 


(4-31) 


where r is the spherical albedo of the atmosphere alone. (Note that the 
relative simplicity of this expression stems from the assumed Lamber- 
tian character of the surface. If the surface were not Lambertian, this 


55 



Introduction to the Theory of Atmospheric Radiative Transfer 

would be a much more complex problem — see Tanre.) The first term 
in the bracket represents flux which has been reflected once from the 
surface, and then scattered back down to the surface. The second term 
represents flux which is reflected upward from the surface, scattered 
down by the atmosphere, reflected back upwards by the surface, and 
finally scattered back downward by the atmosphere. The remaining 
terms are interpreted similarly as multiple reflections and scatterings 
between the surface and the atmosphere. 

The total flux which reaches the surface from the multiple scattering, 
then, is the sum of equations (4-30) and (4-31) 


-frCMo) = rrfioFoTr{-fio) [l + + p^r^ + ■ • -j 

_ nupFoTri-no) 

1 - Psf 


(4-32) 


since < 1. 

Since the surface is assumed to be Lambertian, the total flux 
reflected by the surface is 

^REf(/^o) = PsFrifio) = ^0^ (4-33) 

1 - Psr 

Now, look at the total flux leaving the top of the atmosphere in 
the specific direction cos“^ p. This too is composed of three parts (as 
shown in fig. 4-3): 

1. A component of the incoming solar flux which is directly scattered 
into the direction cos“ ^ p before it reaches the surface 

(4-34) 

2. The total flux received at the surface, reflected by the surface, and 
directly attenuated by the atmosphere 

TrpoFoPaTr{-po) ^_^*/^ ( 4 _ 35 ) 

i^-Psr ^ ^ 


3. The total flux received at the surface, reflected by the surface, and 
diffusely attenuated by the atmosphere 


nflpFoPaTri-po) 
1 - Paf 


t{p) 


(4-36) 
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Figure 4-3. Sketch of some of the ways a given incoming photon can interact with 
the atmosphere and surface and finally escape. 


Thus, writing 

= (4-37) 

we can write the total flux leaving the top of the atmosphere in the 
direction defined by fi as 

^ ^ r,. ^ , ‘^t^oFoPsTr{-IXo)Tr{n) 

Mo) = 7T^o^O^(M» Mo) H ^ (4-38) 

^ ~ Ps' 

Following equation (4-16), we can define the total bidirectional 
reflectance of the atmosphere-surface as 

(4-39) 

^INC 1 “ Ps^ 

(Note again the symmetry r*{p,fio) = r*{fiQ,fj,).) Then, by analogy 
with equation (4-20), if we multiply equation (4-39) by 2p dp and 
integrate over all p, we get the planar albedo of the atmosphere and 
surface system, 

r*(Mo) = r{po) + (4-40) 

where ^ 

fr^2[ pTr{p)dp (4-41) 

Jo 

Finally, the spherical albedo is obtained from equation (4-40) by 
integrating over pq: 

f*^f + (4-42) 

Liou (1980) develops these same relations in a much more rigorous 
way by applying the basic definitions of the R and T functions to the 
RTE. It is felt, however, that the more heuristic approach given here, 
following Tanre, brings the physics of the process more directly into 
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the derivation, and, hence, may be more appealing to the reader, who 
wants to see physically how the various terms react with each other. 

Liou’s development should not be ignored, however, as it permits 
one to derive the same results by a more rigorous manipulation of the 
basic definitions and concepts and, hence, to attain some fluency in the 
use of these more formal statements. In this same context, see also 
section 72 of Chandrasekhar (1960). 

It should be pointed out that for homogeneous atmospheres t{p) — 
but this is not generally true for nonhomogeneous atmospheres; 
i.e., the upward transmission function for a nonhomogeneous atmo- 
sphere is not, in general, equal to the downward transmission. 
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Approximate Analytical Solutions to the RTE 

There are a number of approximate solutions to the various forms 
of the RTE we have developed so far, and considering their simplicity, 
for the Earth’s atmosphere many of them are surprisingly accurate 
when compared with “exact” solutions. The reason for this is that, 
except in the cases of radiation through clouds, heavy fog, or haze, the 
Earth’s atmosphere is optically thin. Many of the approximate solutions 
are based on thin atmospheres, which allow only very low orders of 
scattering to dominate, and thus, when applied to many problems 
or studies in the Earth’s atmosphere, yield numerical solutions which 
compare very favorably in accuracy with much more sophisticated 
“exact” solutions. However, some care must be taken to insure that the 
solutions presented in this chapter are only applied under the conditions 
for which they were derived. Long-term familiarity with, and perhaps 
daily application of such solutions, frequently causes even the expert 
to forget the regions of applicability, so one must be wary of trying to 
apply these approximate results to problems for which the generating 
assumptions are not valid. 

The first two solutions covered in this chapter — the thin-atmosphere 
approximation, and the single- scattering solution— are either applicable 
only to, or are generally more accurate when applied to, a thin atmo- 
sphere; i.e., an atmosphere dominated by low orders of scattering. This 
can occur in an atmosphere of small optical depth, or in an atmosphere 
of large optical depth if its absorption is also large — i.e., ojiy < 1. (See 
the discussion in Irvine, 1968, or Irvine and Lenoble, 1973.) The sundry 
forms of the two-stream solution are applicable to atmospheres of any 
thickness. 

The two-stream solutions, presented later in this chapter, and the 
Eddington solutions of the next chapter, are examples of a frequently 
recurring theme in RTE work, namely, the directional averaging of 
the radiance in order to achieve computationally tractable results. 
The methods of Schuster-Schwartzschild, Sagan-Pollack, and Coakley- 
Chylek all use different directional averaging devices to reduce the 
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radiances in the upper and lower hemispheres to constant parameters 
independent of direction. This process results in a pair of coupled 
linear differential equations with constant coefficients (for homogeneous 
atmospheres), one for the upward intensity and one for the downward 
intensity. 

Thin-Atmosphere Approximation 

The thin-atmosphere approximation is probably the most direct and 
simplest solution of the RTE. (We assume here only the azimuthally 
symmetric case.) It can be obtained directly from the RTE by simply 
assuming that the atmosphere is so thin optically that the derivatives 


dr^, 

can be replaced by their finite difference forms — see equations (5-7) and 
(5-8) and Coakley and Chylek, 1975. 

Start from equation (2-41) 

I I 

with the normalized, azimuthally averaged phase function 

= 1 ( 5 - 2 ) 

Recall from the discussion in chapter 2 that equation (5-1) contains 
both the direct and diffuse radiation components. 

Separate the integral in equation (5-1) into the upward and down- 
ward components 


~yJo V 

dTu I Jo 


60 


Chapter 5 


Now, let us distinguish the upward and downward components of 
intensity by the symbols 


and then write equation (5-3) for each component separately. For the 
upward component 

^ dp' 

dTi/ ^ ^0 

P lUr.,p')P{p,-p') dp' (5-4) 

^ Jo 

To get the downward component, replace p with —p in the first, second, 
and fourth terms of equation (5-4). This is not necessary for ll in 
the third term, as this is the upward component of intensity which is 
scattered downward: 

dlUr.,p) ^ p ll[r,,p')P{-p,p') dp' 

dTu ^ Jo 

P lUru,p')P{-t^, -f^') dp' (5-5) 

^ Jo 

Recognizing that 

p{-p,-p') = p{p,p') 

(see eq. (1-43)) we can write equation (5-5) as 


_^dlh^^ = 4(r.,/i) P ll{r,,p')P{-p,p') dp' 

dtjy I Jo 

P lUr.,t^')Pit^,f^') dp' (5-6) 

2 Jo 

Now, replace the derivatives with the finite difference forms 

/i) ^ ~~ /[/(Oi /^) (fi"7) 

dti/ Ti/ 


and 


dih {ju y M) 


dTi/ 


Ti, 


(5-8) 
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Put equation (5-7) into equation (5-4) and solve for 

4(0, m) = (l - — ) + y— f dfi' 

\ M i fi Jo 

+ -mO V (5-9) 

M ^ JO 

Note the physical significance of the terms in equation (5-9) as sketched 
in figure 5-1. The first term on the right-hand side is the upward 
intensity at in the direction fi which is not scattered— it only 
undergoes absorption along the slant path from Ti^ to Ti/ — 0: 

/ 



Figure 5-1. Sketch illustrating the physical interpretation of the terms of 
equation (5-9). 

The second term is the upward intensity at in the upward 
direction which is scattered into the direction /i, and the third term 
is the downward intensity at in the direction —pf which is scattered 
upward into the direction p. 

Now, we want to eliminate the ll term in the second integral of 
equation (5-9). So we substitute equation (5-8) into equation (5-6) and 

solve for ll{ri^^p) 

4(^1^. m) = 4(0. m) / 4(t'«/.4)'P(-/^.m') dfi' 

\ M / M ^ ^0 

+ — dti' (5-10) 

H 2 Jo 

The terms in equation (5-10) have a similar interpretation to those of 
equation (5-9). If equation (5-10) is substituted into equation (5-9) and 
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only first-order terms in Ti^ are retained, then 
and equation (5-9) becomes 

;'{0,m) = f 1 - f' V 

\ M / p Z JQ 

+ -V V) V (5-11) 

M ^ Jo 

We can get the reflection and transmission coefficients of a thin 
atmosphere directly from equation (5-11) (see Coakley and Chylek, 
1975). To get the reflection function, assume a solar beam incident 
on the top of the atmosphere 

27t/^(0,)[/o) = 7rFo<5(/i - no) (5-12) 

(the factor 27 t arises from the azimuthal integration) and assume the 
incident diffuse radiation at the bottom of the atmosphere 

27r/2(rj/,/i) = 0 (5-13) 

Then equation (5-11) becomes for this case 

^/I(o./^) = 

The reflected flux is given by 

rl ^ /* 1 

Fl{n) = 27t / /i/^(0,/z) dn = -TyCbuFo / P{fi,~no) dn (5-15) 

Jo 2 Jo 

and from equation (4-17) the planetary albedo is 

^(mo) / P{F^-Fo)dn (5-16) 

2 no Jo 

We can use equation (5-11) to find the transmission function, T{no), 
even though equation (5-11) describes the upward intensity component, 
by the simple artifice of letting the solar flux impinge on the bottom of 
the layer and examining the flux emerging from the top of the layer. 
Note that we use T{no) rather than t{no) to denote the transmission 
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function, as here we are determining the total transmission function and 
not just the diffuse component. 

We have the boundary conditions 

= t<-FoS{h - fio) (5-17) 

and 

27t/^(0, //) = 0 (5-18) 

Put these into equation (5-11) to get 

4(0, M) = (i - 6{fj. - w,) + 

Again the flux is given by equation (5-15) 

fJ (0) = 27T/io (l - ^) y + y P{fi\ Mo) dfi’ 

The transmission function is gotten from equation (4-19) 

T{^^) = (l - ^) + r P(4,//o) V (5-19) 

\ M/ M 2 Jo 

The first term on the right-hand side of equation (5-19) represents the 
contribution of the direct transmission 



and the second term, which is analogous to equation (5-16), is the 
transmission function for the diffuse term. 

We can use the normalized property of the phase function, equa- 
tion (1-27), to write this in another form. From equation (1-27) and 
azimuthal symmetry, we have 

= 1 

or 

^ 1° P{^i, 4) d^l' + F(/i, 4) = 1 
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and since 

P{fi,-fi') = P{-^i,li') 

we can write the integral form of equation (5-19) as 

^ ^ Pin', no) = ~f^o) dn' 

and hence write equation (5-19) as 

T{no) = 1 - — [i - wz. + ^ Pin, -w) dn] (5-20) 

MO L ^ Jo 

Equations (5-16) and (5-20) show that for a thin atmosphere, both 
the albedo and the transmission functions are linear functions of optical 
depth. 

Define ^ 

= ^ / P{p,-m)dn ( 5 - 21 ) 

2 Jo 

Then from equation (5-16) 

rino) = —Cj^0ino) (5-22) 

Mo 

And from equation (5-20) 

^(mo) “1 ^ [1 — a)i/ -f oJi/f3{po)] (5-23) 

Mo 

Finally, if we define 

= / /3(mo) ^MO (5-24) 

Jo 

we can get the spherical albedo and spherical transmission function by 
using equations (5-22) and (5-23) in equations (4-24) and (4-25) 

f = (5-25) 

f - 1 - 2Ti. (1 ~ oJi. -h (5-26) 

The quantities /3 (mo) ^ ^sed quite extensively in the litera- 

ture, especially that pertaining to the derivation of approximate solu- 
tions to the RTE. The quantity /3 (mo) is the bachcatter fraction for a 
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beam of radiation entering the atmosphere at the angle cos“^ pQ. This 
is geometrically proportional to the fraction of the total surface area 
of the phase function above the horizontal plane through the scatter- 
ing center, as sketched in figure 5-2. The quantity ^ is the integrated 
backscatter fraction over the whole range of entry angles. Azimuthal 
symmetry has been assumed throughout this section, and hence, also 
in the definition of /3{fio) 



Figure 5-2. Interpretation of the backscatter fraction. The phase function is the 

Henyey-Greenstein type for small asymmetry parameter, g. 

Wiscombe and Grams (1976) discuss these backscatter fractions in 
detail and give integral methods of evaluating them for general phase 
functions. Table 5-1 shows computed by their method for various 

values of g (asymmetrical parameter) and fiQ^ and table 5-2 shows values 
of Both tables are for the Henyey-Greenstein phase function. The 
table data are also plotted in two accompanying figures (figs. 5-3 and 
5-4). 

The tables reflect one’s intuition about the behavior of /?(//q) 

For isotropic scattering, /?(//q) = /? = 1/2 for all fiQ (one-half of the 
radiation is scattered forward and one- half is scattered backward for any 
entry angle). For very elongated phase functions [g near 1), most of the 
radiation is scattered in the forward direction. Hence, both /?(//) and 
approach zero (very little backscattering). For very low incidence angles 
near 90 degrees (/.f 0) a somewhat higher fraction is backscattered 

than for near-normal incidence angles, j3{pi = 1) > /3(/Lio = 0). 

The reflection and transmission functions from equations (5-22) and 
(5-23) are compared with some exact computations using the doubling 
method (Liou, 1973) in figures 5-5 and 5-6. Note that, as expected, 
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TABLE 5-1. VALUES OF /?(mo) vs. g FOR VARIOUS po 
[For the Henyey-Greenstein phase function] 


g 

o 

II 

o 

d 

II 

O 

fiQ - 0.3 

o 

11 

o 

o 

II 

p 

d 

0.00 

0.500 

0.500 

0.500 

0.500 

0.500 

.05 

.496 

.492 

.489 

.485 

.481 

.10 

.492 

.485 

.477 

.470 

.462 

.15 

.489 

.477 

.466 

.454 

.443 

.20 

.484 

.469 

.454 

.438 

.423 

.25 

.480 

.460 

.441 

.422 

.403 

.30 

.476 

.451 

.428 

.404 

.382 

.35 

.471 

.442 

.413 

.386 

.360 

.40 

.465 

.431 

.398 

.367 

.337 

.45 

.459 

.419 

.381 

.346 

.313 

.50 

.452 

.405 

.362 

.323 

.288 

.55 

.443 

.390 

.341 

.298 

.262 

.60 

.434 

.372 

.318 

.272 

.234 

.65 

.421 

.350 

.291 

.243 

.205 

.70 

.406 

.324 

.260 

.211 

.174 

.75 

.385 

.292 

.225 

.177 

.144 

.80 

.356 

.252 

.185 

.142 

,113 

.85 

.314 

.202 

.141 

.105 

.085 

.90 

.247 

.141 

.094 

.069 

.053 

.95 

.141 

.072 

.046 

.033 

.026 

1.00 

.0 

.0 

.0 

.0 

.0 


the solutions (5-22) and (5-23) for the thin atmosphere (r = 0.0625) 
show better agreement with the exact calculations than do those for 
the thicker atmosphere (t = 0.25). For both thicknesses, the agreement 
is also better for steep incidence angles (/tq 1) than for shallow 
incidence angles (/io ~ 0)? because for steep entries there is less chance 
for multiple scattering to occur. 

Single-Scatter Solution 

The single-scattering solution to the intensity equation is probably 
the next simplest solution to the RTE. In this solution, we permit the 
incoming solar radiation to be scattered only once, and compute the 
resulting upward and downward intensities resulting from this single 
scatter. 

Many phenomena involving atmospheric scattering can be ade- 
quately represented by the single-scattering approximation, the most 


67 




Introduction to the Theory of Atmospheric Radiative Transfer 

TABLE 5 - 1 . Concluded 


g 

HQ = 0.6 

HO - 0.7 

HO = 0.8 

HO = 0.9 

PO ~ 10 

0.00 

0.500 

0.500 

0.500 

0.500 

0.500 

.05 

All 

.474 

.470 

.466 

.463 

.10 

.455 

.440 

.433 

.425 

.423 

.15 

.432 

.421 

.410 

.399 

.389 

.20 

.409 

.394 

.380 

.367 

.353 

,25 

.385 

.368 

.351 

,334 

,319 

.30 

.361 

.340 

.321 

.303 

.286 

.35 

.336 

.313 

.292 

.273 

.255 

,40 

.311 

.286 

.263 

.243 

.225 

.45 

.284 

.258 

.236 

.215 

.197 

.50 

,258 

.231 

.208 

.188 

.171 

.55 

.230 

.204 

.181 

,163 

.147 

.60 

.203 

.177 

.156 

.139 

.124 

.65 

.175 

.151 

.132 

.116 

.103 

.70 

.147 

.125 

.109 

.095 

.084 

.75 

,119 

.101 

.087 

.076 

.067 

.80 

.093 

.078 

.066 

.058 

.051 

.85 

,067 

.056 

.048 

.041 

.036 

.90 

.043 

.036 

.030 

,026 

.023 

.95 

.021 

.017 

.014 

.012 

.011 

1.00 

.0 

.0 

.0 

.0 

.0 


notable exceptions being the scattering characteristics of clouds, heavy 
haze, and fog, and possibly heavy aerosol concentrations. The extinc- 
tion coefficient for background aerosols in the stratosphere, for example, 
is of the order of 2 x 10^"^ km. Thus, the mean free path for stratospheric 
aerosol extinction is of the order of 5000 km, and the single-scattering 
solution should suffice for all but the most shallow solar flux entry an- 
gles, along which the possibility of more than one scatter might take 
place (see Buglia, 1982). In the troposphere, a dear-day extinction 
coefficient might be of the order of 2 x 10”^ km, giving a mean free 
path of the order of 50 km, so that even here the single-scatter solution 
might be used for some problems. In a heavy fog or haze, the extinction 
coefficient might be of the order of 1 to 10 km, and obviously one could 
not try to use the single-scatter solution under these conditions. 

We start with the formal solutions for the upward and downward 
intensities, equations (3-5) and (3-7), which we now write, dropping the 
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Chapter 5 


9 

p 

9 

0 

0.00 

0.500 

0.55 

0.283 

.05 

.481 

.60 

.261 

.10 

.462 

,65 

.238 

.15 

.444 

.70 

.214 

.20 

.425 

.75 

.188 

.25 

.405 

.80 

.161 

.30 

.386 

.85 

.131 

.35 

.366 

.90 

.098 

.40 

.346 

.95 

.058 

.45 

.326 

1.00 

.000 

.50 

.305 




V subscript, as 

I(r,^,,4>) = (5-27) 

+ f J{t ,-n,<j>)e~^'^ (^28) 

JQ ^ 

The source function, J(r, is the singly scattered incoming solar 
radiance, which we write as 


J(r,/x,0) ttFoc 


An 


and which is the product of three terms: 


(5-29) 


ttFoc-^/^o 




the incoming direct solar intensity attenu- 
ated to the level r 

the single-scattering albedo; i.e., the frac- 
tion of the incoming solar radiance which 
undergoes scattering 

the fraction of the scattered radiance which 
is scattered from the direction (—po^^o) 
into the direction (p, 0). 
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0 .2 .4 ,6 .8 1.0 

^0 


Figure 5-5. Comparison of the thin-atmosphere approximation, equations (5-22) 
and (5-23), with the exact (doubling) method for two values of the optical 
depth, r = 0.0625 and r = 0.25. The Henyey-Greenstein phase function was 
used with g — 0.75, cDq = f-0- 
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We assume as boundary conditions 


= 0 (5-30) 

= Q (5-31) 


i.e., no diffuse radiation enters the top or bottom of the atmosphere. 
From equation (5-27), with equation (5-29) and the boundary 
conditions 




dr' 


(5-32) 


In particular, at the top of the atmosphere r = 0, and we get 




4{fi + fio) 





(5-33a) 


Comparison of equation (5-33a) with equation (4-8) shows that we can 
write the reflection coefficient for single scattering as 


4 fl+flQ 



(5-33b) 


In a similar way, we get the downward component of the intensity 
by substituting equation (5-29) and the boundary conditions into 
equation (5-28) 


■kFqP{-h, 4> : -fio,<l>o) f e e 

Jo 


dr' 

(J- 


uj 

= -tPoP{-(^, (f> ■ -fM), <l>o) 

4 fi 


Here, we must distinguish between two cases: 

1. /z = /zq, and 

2. fiQ. 
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For case 1 we get immediately 

h{T,~n,4>) = (5-34) 

and for case 2, 

/ 2 (r, -M, 4>) = ^ 0 : -MO, 0o) (5-35) 

Emerging from the bottom of the atmosphere 

7i(r*,-M, 0) = ^PoP(-M>^ : -MO,<^o)— (5-36a) 

= 7 ■ -Mo,0o) [e”’’ -e“'^ 

^ (5-36b) 

Comparing equation (5-36b) with equation (4-9) gives the diffuse trans- 
mission coefficient for single scattering 

< 2 (m,m) - (a^ W)) (5-37a) 

The direct component of the transmission coefficient is, of course, 

g-r'/MO 

For case 1 (m = /^o), the diffuse part of the transmission function is 

ti(M,Mo) = = ^o) (5-37b) 

't Mo 

We can now easily show that for a thin atmosphere, t* •<C 1, 
equations (5-32) and (5-35) reduce to the thin-atmosphere solution 
derived earlier. For r* < 1, we get from equation (5-32) 

/(r,M,0) = 7 ^^^P(m,*^ : -MO,<Ao) (5-38) 

4 M 

Assuming azimuthal symmetry the upward flux at r = 0 is 

F^(0) = 27t / mt^O— P(m, -Mo) ^^M = TrCjFQT* l3{fio) 

^ 4 M 
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which is identical to equation (5-15), derived from the thin- atmosphere 
assumption. 

The direct component of the downward flux is 

ttFo/zq ~ (5-39) 

The diffuse downward component of intensity becomes, from equa- 
tion (5-36b) 


7(r*, j : -/^0^</>o) U 1 + — ] = -fiQ) 

I pt MoJ 4 ^ 

and the diffuse component of the downward flux becomes 

F^{t*) = 2n [ f^jFc — P{-fi, -no) dfi = [1 - j3{no)] 

Jo F 

The total downward flux is thus 

F^{t*) = TrFono fl - — ) + 7rFowr*[l - p{no)] 

\ FoJ 

From this, the total transmission function emerges as 
T{no) = 1 - — [1 - d) + Cb^ino)] 

FO 

which is identical to equation (5-23). Note the difference between the 
thin-atmosphere and single-scattering solutions. The single-scattering 
solution makes no assumptions about the thickness of the atmosphere — 
it only assumes that the photons are scattered only once. 

Two-Stream Solutions 

The two-stream solutions are, in general, arrived at by writing the 
RTE for the upward and downward components separately, and assum- 
ing that the upward intensity is constant over the upper hemisphere 
and independent of the angle //, and that the downward component 
is a different constant over the lower hemisphere, also independent of 
fi. The differential equations each involve an integral, and the method 
of approximating the integral leads to a set of two linear differential 
equations. These equations have constant coefficients when applied to 
homogeneous atmospheres, as they are here. 
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We consider in this chapter three such methods of approximat- 
ing the integral mentioned above. These lead to well-known ap- 
proximate solutions — the Schuster-Schwartzschild approximation, the 
Sagan-Pollack approximation, and the Coakley-Chylek approximation. 
This third form, the C-C approximation, appears to have a slightly 
less rigorous formulation than the others, but it retains the dependence 
of the solution on the solar incidence angle, /io- When appropriately 
applied, these equations all give numerical results which are in good 
agreement with other solutions. 

The differential equations resulting from these three approaches are 
identical in algebraic form, the only differences appearing in the terms 
making up the constant numerical coefficients. 

These three sets of two-stream approximations will be derived in 
this subsection. The reader is invited to examine the excellent review 
article by Meador and Weaver (1980), in which a number of well-known 
approximations, including other than the classical two- stream solutions 
developed here, are discussed and compared. Meador and Weaver 
neatly identify the theoretical thread common to all of these methods 
and show that they all reduce to the same algebraic form, except for 
the grouping and definition of some constant algebraic parameters. The 
paper by Lyzenga (1973) is also worthy reading, in that he shows that 
the Sagan and Pollack formulation can in fact be rigorously derived by 
assuming a two-point Gaussian quadrature formula to approximate the 
integrals mentioned above. He also shows that a single transformation 
relates the two-stream and the Eddington approximation discussed in 
the next chapter. Lyzenga’s approach is used below to derive the Sagan- 
Pollack equations. 

The two-stream analysis is applied to the azimuthally symmetric 
form of the RTE, for either the total intensity (direct plus diffuse, 
eq. (2-41)), or for the diffuse component only (eq. (2-50)). We will 
not derive in detail all the possible combinations here, as the repetition 
would serve no purpose, but will derive one total intensity solution and 
one diffuse intensity solution. Some limited numerical comparisons will 
also be given. 

The Schuster-Schwartzschild (S-S), the Sagan-Pollack (S-P), and 
the Coakley-Chylek (C-C) equations can all be reduced to the same 
differential equation form. These three forms will be derived separately 
below, and the general form of the solution given. 


Schuster-Schwartzschild (ref, e.g,^ Ozisik, 1973). Start with 
the form of equation (2-41) and write the upward and downward 
components separately, as in equations (5-4) and (5-6). We will drop 
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the subscript u from the /, the r, and the w, but must keep in mind 
that these developments are for monochromatic radiation only. 


-/i- 


d/^(T,/x) 


dr 




f (5-41) 

Multiply equation (5-40) by dp and integrate from 0 to 1 

J j dn j f[r,ti)P[n,n')dn' 

f (^ 42 ) 

Interchange the order of integration in the last two terms 

fi) dn - J i\t, iji') dfi J P{fi,fi')dn 
l\r,ti')P(n,-n')dn' (5-43) 

Since by equation (1-33) and the symmetry in Pyp\ 


we have with the definition of equation (5-21) 

f /^{r,M')V (5-44) 

Jo 

To this point, equation (5-44) is exact — at least insofar as the differ- 
ential equation (5-40) is exact. Now, we make one of the approximations 
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mentioned above — we assume I to be independent of p in each hemi- 
sphere, so that /^(r, /z) — ► /^(r) (and similarly for /J^). This gives the 
S-S approximation 

dfi « (5-45) 

and we write equation (5-44) in the form 

1 = I\t)-oj[I - ^(//)]/^(r) - w/?(/x)/^(r) (5-46) 

Z ar 

In a similar manner we can develop equation (5-41) to the form 


(5-47) 

Equations (5-46) and (5-47) are the S-S form of the differential equa- 
tions describing the upward and downward components of the total 
intensity fields in the two-stream approximation. 

Sagan-Pollack (Sagan and Pollack^ 1967; see also Lyzenga, 1973). 

Again we start with equation (2-41) with the subscript u dropped from 
/, r, and u 

- 1 1 ', 

Lyzenga argues that instead of taking and to be some average 
value of I over their respective hemispheres, it is more appropriate to 
be guided by the two-point Gaussian quadrature method of numerical 
integration, and take for the average value of / that value which 
would be obtained if we solved equation (5-48) along the ray given by 
pL = ±l\/3. This conclusion can be substantiated more formally from 
equation (2-41), if we replace the integral with a two-point Gaussian 
quadrature and write our equation for each ray separately. This gives 
the pair of equations 

= ( 5 - 49 ) 
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and 


- 1 /_‘ 

Apply the Gaussian two-point quadrature formula to the integrals in 
equations (5-49) and (5-50) 


(5-52) 

If we use the two-term expansion of the phase function, equation (2-42), 
we get 

and hence, the integrals in equations (5-51) and (5-52) become, 
respectively, 

and 

and equations (5-49) and (5-50) become 


^ (r) - w( 1 - 6)/T (r) - Cjbl^ (r) 


(5-53) 


and 


in both of which 

6=i(l-|) (5.55) 

These are the S-P forms of the two-stream equations. We can further 
evaluate u)i in terms of the more familiar aisymmetry factor, which 
is the first moment of the pheise function (see eqs. (1-46) and (1-47)) 


g 


1 

2 




dfi 


(5-56) 
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and with equation (2-31) we get 

N 


1 

1 ^ /-I 

= oYl^j 

S=i 

\ ^ rl 

= ol2^J Piif^)PjMd(x = 0 
S=i J-i 

= 2"‘/_/ ‘'"“y 

Hence, we can define b in the more familiar way 

t = 5(i-») 


(J = l) 


(5-57) 


(5-58) 


Coakley-Chylek (Coakley and Chylek, 1975). The C-C form of the 
two-stream equations is found directly from equations (5-4) and (5-6) 
simply by assuming that and 7^ are independent of n, and using 
the definition of equation (5-21). This gives immediately the pair of 
equations 

- 7^ (r) - (D[l - /?(/x)]7t (r) - (5-59) 

- w/3(/i)7^(T) -a)[l -P{ix)]I^{t) (5-60) 

Solution of the Two-Stream Equations 

Comparison of equation (5-46) with equation (5-47), equation (5-53) 
with equation (5-54), and equation (5-59) with equation (5-60) shows 
that they all can be put into the same algebraic form 

^1^7^ =7^(r)-w(l-7)7^(r)-w77^(r) (5-61) 
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= I^{t) ~Cj{1 - 7)/^(r) (5-62) 

in which we have the correspondences recorded in table 5-3. 

TABLE 5-3. COMPARISON BETWEEN AND 7 
FOR THE THREE SOLUTION METHODS 


Solution 

method 


7 

s-s 

1/2 


S-P 

l/>/3 

b 

c-c 




It is comforting to note that equations (5-61) and (5-62) satisfy our 
physical intuition. For example, let us look at equation (5-61) at some 
altitude. If we increase z to z-\- dz^ then r decreases by dr. If we write 
equation (5-61) in the form 

= (1 ~ o))/^(r) -h 0)7/^ (r) - 
dr 

we see that I^{z -h dz) is reduced by the first right-hand-side term, the 
absorption of the upward radiance between ^ and ^ -h dz^ as well as by 
the second right-hand-side term, the radiance backscattered out of the 
upward beam, and that it is increased by the last term, the part of the 
downward beam which is backscattered in the upward direction. 

Equations (5-61) and (5-62) can be solved by any number of stan- 
dard techniques. We select here the operator approach. Put equa- 
tions (5-61) and (5-62) into the form 

- [1 -w(l -7))| (5-63) 

+ [1 - (^(1 - 7)]| (5-64) 

Solve equation (5-64) for /^(r). Substitute into equation (5-63) and 
expand the operator to get 
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where 

e = a^-p^ 

(5-66) 

1 - w(l -7) 

a = 

Ml 

(5-67) 

II 

(5-68) 

Equation (5-65) then solves immediately as 


/T (r) = Ae^^ + Be~^^ 

(5-69) 

where A and B are constants of integration. Put equation (5-69) into 
equation (5-63) and solve for 

I^{t) = Awe^^ + Bve~^^ 

(5-70) 

where 

a - ^ 

(5-71) 

and 

o-t-^ 

P 

(5-72) 

Finally, apply the boundary conditions 


11 

(5-73) 

/^(r*) =0 

(5-74) 


where t* is the total optical thickness of the atmosphere. Solve the 

resulting expressions for A and B to give the intensity solutions in final 

form _ .... , * ,■ 

g-C(r - t ) _ -r) 


lHr) = Io 


/^(r) = /o 


we 


-€(r‘-r)_^ee(r*-r) 


(5-75) 


(5-76) 


we — ve^’’^* 

From these, the reflection and transmission functions are found to be 


I^0) _ 

Iq — ve^^* 


(5-77) 
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and 


X = = w- V 

/q 


(5-78) 


Thus, if we for instance use the S-P parameters from table 5-3, 
the resulting expressions may be algebraically different from those of 
Sagan and Pollack (1967), but the results will be numerically identical 
to theirs. 


Solution for conservative scattering, a; = L The case of conservative 
scattering, d; = 1, cannot be found directly from the solutions of 
equations (5-75) and (5-76) simply by setting cD = 1, because ^ = 0 
for u) = 1, and the solution falls apart. The conservative scattering case 
must be gotten by starting with the differential equations (5-61) and 
(5-62) and setting d; = 1 

= llHr) - iI^{t) (5-79) 

-Ml = 7^^’-) - itHr) (5-80) 

The right-hand sides of these equations are identical and must, there- 
fore, be constant 

7 /^ (r) - 7 /i ( 7 -) = M (5-81) 

M is a constant. Equation (5-79) has the immediate solution 

iHt) = — + B (5-82) 

/^1 

where B is again a constant of integration (not the same B as we just 
used earlier). Similarly, equation (5-80) gives 

iHt) = — + B' (5-83) 

/^1 

Substitute equations (5-82) and (5-83) into equation (5-81) to evaluate 
the constant M. Since M is constant for all r, it is most conveniently 
evaluated at r = 0. This gives 


M = B') 
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and equations (5-82) and (5-83) are rewritten as 

/T (r) = B f 1 + — ) - —B' (5-84) 

V fJ'iJ 

(r) = +(i-^)b' (5-85) 

The constants B and S' are evaluated from the boundary conditions 
of equations (5-73) and (5-74), and hence, the final solutions for the 
two-stream conservative scattering case can be written as 


/^(r) = Iq 


— zpt — 


(5-86) 




1--^ 
l + lll 


The reflection and transmission functions become 


(5-87) 




T = 


/t(0) 

111 

Ml 

(5-88) 

7o 

^ Ml 


1 

(5-89) 

Jo 

l + Tl 


Note in this case that R + T = 1, as there is no absorption in the case 
of conservative scattering. Also, in the limit as r* — > oo for d) = 1, 


R{t* — > OO) = 1 (d) = 1) 

T{t* oo) = 0 (w = 1) 

Since there is no absorption, all of the incoming radiation eventually 
escapes from the atmosphere [R = 1), and obviously, nothing passes 
through the infinite optical depth (T = 0). 

The two-stream flux equations are simply obtained from the radi- 
ance solutions 


F^(r) = 2n f /r/^(r)dp = 27T/ii/^(r) (5-90) 

Jo 
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F^{r) = 27t / ptl^{r)dp = 2nfiil^{r) (5-91) 

^0 


Solution for diffuse component only. The two-stream solutions 
presented thus far have involved the total intensity — the diffuse plus 
the direct components. Liou (1980) presents a solution for the diffuse 
component only. The solution closely parallels the development given 
above, so just the barest outline will be given here. Start with equa- 
tion (2-51) and evaluate the integral by the Gauss method: 


r\ n 

J ^ f{x) dx^ (5-92) 


J=-n 


(e.g., Abramowitz and Stegun, 1970), where the weights 


^ 1 P2n{x) 


dx 


and the Xj are the zeros of the even-numbered Legendre polynomials. 
Also, we have 


n 



j=-n 

and thus can write equation (2-51) for a ray defined by Pi as 


l—O j= — n 

- (5-93) 

in which we have used the property of Legendre polynomials 

Pni-X) = i-irPnix) 
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For the two-stream solution, we take i = ±1. Then since pi = 
3“^/^ and a\ — a-i (Abramowitz and Stegun, 1970) we can break 
equation (5-93) into two equations 


=/^(r)-w(l-6)/T(r)-w6/^(r) 

dr 


(5-94) 


'Ml =/^(r) — u;(l — 6)/^ (r) — Cbhl^ (r) 

dr 

- + 3gfiofn)e~^/f^ 


(5-95) 


Proceeding as before, writing equations (5-94) and (5-95) in operator 
notation, we find after some messy but straightforward algebra that 




where the following definitions hold: 


u = 


1 — a 


V = 


1 -}- n 


A = 


a — /? 



(5-96) 

-1- Ae-^/^o 

(5-97) 

1 ._« + /? 

„ 2 _ 1-w 


€ — 


1 “ Qg 


a 


1 - fc2^2 


0 = 


Z2fJ-l 

1 - A;2/x2 


2 _ (1 - <:^g)(i - <^) 

/v 


-(1-Wg)(s -t-S+) . (s -5+) 

Zl = 2 1 

Ml MMl 

— (1 — u;)(5'~ — s~^) ( 5 “ -(- 5 ^) 

2^2 = — S h 

Pi POPl 

= ^Fo(l ±3gM0Ml) 

The application of the boundary conditions usually used with the diffuse 
component 

/i(0) = /T(t*) =0 
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gives the constants 


eve 

y2 ^— kr * _ y2 ^ kr * 


Xve^'^ — eue 

U^e~kr* _ y2^kr* 

With the fluxes determined from equations (5-90) and (5-91), the 
planetary albedo is given by equation (4-16) with the incoming col- 
limated flux by equation (4-14). From these, after some algebraic 
reduction, 




ttFo/xq 


+ uv 


(2!-^) (.'*'• (5-98) 


in which 


Gi 


ojpo 




1 - fc2^2 

^ 1 - fc 2^2 [sir + -<^)] 


D = 


The diffuse transmission coefficient can similarly be written as 


i(/xo) = 


t^FqUo 


2/^1 / ( G2 — Gi^ 

= j 




- 2/-1 

The total transmission coefficient, diffuse plus direct is, of course, 

T’(mo) = t{^iQ) + e-^‘/M0 (5.99) 
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If we consider the special case of an infinitely thick atmosphere, 
r* -4 00 , we find that in the limit, the albedo of an infinitely thick 
atmosphere becomes 


r(/zo) = 2/ii 



(5-100) 


Note that this is not unity, as we found before for the conservative 
atmosphere case. Here, Cj ^ I and there is absorption present, so not 
all of the incident radiation escapes from the top of the atmosphere. 

In the next chapter, results from equations (5-98) to (5-100) will be 
compared with those from the Eddington approximation. 

The conservative solution, in which a) = 1, also proceeds from 
equations (5-94) and (5-95) 




dl^ (r) 
dr 


6/^(r) — bl^{r) — s e 


(5-101) 


/^l 


dl^{r) 

dr 




the solutions to which are 


(5-102) 




/i(r) = Si 



_ br 
- B 2 — 
Ml 

- 1 - 

(5-103) 

hr „ / 
— +^2 
Ml ' 

< Ml/ 

+ l3iFoe-^/f^ 

(5-104) 


in which 


Bi = 


1-t- 


br* 


B 2 = —0iFq 


ai — (l - SgfiQfii - 26 —) 

= - j— (1 +3gfiQfj,i - 1 - 26 —) 

4^1 V 
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The fluxes, as before, are given by equations (5-90) and (5-91). The 
upward flux at the top of the atmosphere is 


f\0) = 


2nfioniFo 




1 + 


6r* 

/^1 


(5-105) 


and hence, the planetary albedo becomes 

^ ^ (l - 3gMopi-26^) (l-e-^V^o) 


^(mo) = 


1-h^ 

/^1 


(5-106) 


As it should, this also approaches 1 as r* becomes infinite. 
Similarly, the downward flux at the surface, r = r*, becomes 


F^{r*) — 27r/xi 


B1— + B2 

Ml 


(1 - ^) -h 


and since the scattering is conservative 

i(Mo) = 1 - r(/Lio) 


(5-107) 


from which the diffuse transmission function can be found. 

There is a very interesting extension of the Schuster two-stream 
method to n streams in a paper by Acquista, House, and Jafolla (1981). 
No azimuthal symmetry is assumed, and instead of just considering an 
upper and a lower hemisphere, as done above, the authors break each 
hemisphere into an arbitrary number of nonoverlapping patches and 
compute the radiance stream for each patch. No numerical data are 
presented, but it is reasonably claimed that computational costs are 
substantially less than for other more nearly exact methods for a given 
accuracy, and great flexibility in the choice of patches is claimed, for 
minimizing computational costs for particular applications. 
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The Eddington Solution 

The basic differential equations describing the Eddington approxi- 
mation are derived somewhat differently from those of the two-stream 
solutions. However, the mathematical technique leading to the solution 
is quite similar, and hence, will not be given in the same detail as in 
the preceding chapter. 

The Eddington solution begins by assuming that the intensity, 
instead of being approximated by different upward and downward 
constants, can be approximated by a linear function of of the form 

m) = (6-1) 

where Iq and /i are functions of r only, and not of //. Now, if we apply 
the same two-stream approximation to equation (2-51), i.e., let N = 1, 
we get 

-/(r,/x) - I [c:>oPo(m) l'^PQ{fi')nr,Adti' 

+ ^iFi(m) d(i'] 

- [CjoPo{fi)Poi-tio) + ojiPi{fi)Pi{-fto)] 

4 

or 

Y [/ dn' + Cjiti dfi' 

( 6 - 2 ) 

which is identical to assuming a two-term expansion of the phase 
function directly in equation (2-50) 

P{n,^i') « l+wi/x/x' 
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The Eddington method for solving the RTE also leads to a two- 
stream type of solution, but as indicated earlier, and discussed in detail 
below, begins in a completely different manner. It will be seen that the 
intensity boundary conditions cannot be completely satisfied at either 
surface, although the flux boundary conditions can be at least formally, 
if not exactly, physically satisfied. For this reason, the Eddington 
solution is more accurate for very thick optical depths. It also uses a 
two-term expansion for the phase function, resulting in a solution that 
is most accurate for scattering, which is close to isotropic. This occurs 
well inside an optically thick atmosphere, after multiple scatterings have 
taken place (see Irvine, 1968). Multiple scattering deep in the interior 
effectively smooths out the phase function in that the sharp maxima 
and minima usually present in the phase function disappear, and the 
scattering becomes more nearly isotropic. 

The integrals in equation (6-2) can be evaluated with the use of the 
Eddington assumption, equation (6-1) 

/(r,/z')V = 2/o(t) 

= ^hir) 

and with equation (5-57), equation (6-2) can be written as 

[^(7-) + W] + f^h{r) - w [/o(r) + gfiliir)] 

(6-3) 

Equation (6-3) can be broken up into two equations for /i(r) and 
Iq{t) as follows: multiply equation (6-3) by dp and integrate from —1 
to ”f“l. 

= 3(1 - u)Iq{t) - (6-4) 

Now, multiply equation (6-3) by pdp and integrate over the same limits 

= (1 - w</)/i(r) -I- ^u;Fog/jioe~''^'^ (6-5) 

We now again have two coupled linear ordinary differential equations 
with constant coefficients. These are solved the same way as for the 
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two-stream solution, and give 

/ o ( t ) = + Be~^^ + 

/i(r) = aAe^'^ - aBe~'^^ + 
in which the following definitions hold 


^2 _ 3(1 — d)) 
I -ug 



= 3(1 — w)(l — Cjg) 


and 


Go = 1 — ^g + g 
Gi = 3(1 —uj)gnQ H 

Note that here we cannot apply the boundary conditions 


( 6 - 6 ) 

(6-7) 


/i(0)-0,/^(r*)-0 

since to do so would result in two equations in four unknowns. The 
reason for this is that equation (6-1) is really the first two terms of a 
Legendre polynomial expansion for the radiance, /(r, /x) 


7(r,/x) = /o(^)^o(m) + h{r)Pi{fi) 

and we cannot depict a constant function (constant zero flux at all entry 
angles ^o) with only a two-term expansion. However, we can at least 
formally apply the boundary conditions to the flux form of the solution. 
Using equation (6-1) in the flux definition 


r 2 ' 

F^(r) = 27ry^ ^ v Io{t) + -Ii{t) 

fO r 2 ' 

F^(r) = 27ry n)dfi = ir Iq{t) - 

Insert equations (6-6) and (6-7) into equations (6-8) and (6-9) 

(r) = -I- Bue~^^ -I- 


(6-8) 

(6-9) 


( 6 - 10 ) 
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F^{t) ^ Aue’^^ + Bve~'^^ + (6-11) 

where 



Then, from the flux boundary conditions 

F^t*) = 0 F^(0) = 0 


vee ’’ /^o — 'yue 

u^e~kT* _ y2^kT* 
vfie^kr* _ y2^kr* 

Just £is we did in chapter 5 for the two-stream case, we can write 
the reflection and diffuse transmission coefficients for the Eddington 
solution as 

?*(mo) = - L2UV ^ j + Li (6-12) 


e(/io) = 
where 


— ~ L2 (i/^ — u^) + Liuv ^ (^0 — e ^ (mo ^ j ^ /#^0 

(6^13) 


(jj 


Li =- 


3(1 - oj)g(4 -I - 1 - §(i - + g)Mo 

1 - k'^nl 


T 

i2- 2 


3(1 — oj)gfiQ -t- 1 -l- §(1 — Cjg + g)no 
1 - k'^nl 


D 




2 icT* 
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For a semi-infinite atmosphere 

,(«) = (^ 14 ) 

The case of conservative scattering, d; = 1, must again be handled 
separately. Start with equations (6-4) and (6-5) and set d) = 1 


(6-15) 

dr 4 

and 

= (1 - g)h{r) + (6-16) 

Equation (6-15) immediately integrates to 

h{r) = ^Fotioe-^^>^^ +K (6-17) 

where if is a constant of integration. Substituting equation (6-17) into 
equation (6-16), and integrating, yields 

Hr) = + (1 _ g)KT + H (6-18) 


where H is another constant of integration. 

If now equations (6-17) and (6-18) are put into the flux equa- 
tions (6-8) and (6-9), and the same boundary conditions are applied, 
we get finally 


H = 


:K 


(6-19) 


with 


K = 


[l + Imp -h (i - Imp) 


3(l-g)r* -t-4 


( 6 - 20 ) 


From these results, the albedo for conservative scattering in the 
Eddington approximation becomes 


where 


^•(mo) = 1 - 


2L (r*,/io) 
3(1 -g)r* +4 


L{t*,ho) = H- + (l - ^Mo) e 


( 6 - 21 ) 

( 6 - 22 ) 
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Two-Stream and Eddington Solutions for Semi-Infinite 

Atmospheres 

For an infinitely thick atmosphere, r* — ► oo, the two-stream albedo 
is given by equation (5-100) and the Eddington solution by equa- 
tion (6-14). 

The two-stream and Eddington solutions for a semi-infinite at- 
mosphere are compared in figure 6-1. (The figure also shows other 
solutions — including the delta-Eddington which is discussed below.) 
The particular numerical values were selected to permit comparison 
with numerical results in Irvine (1968). We can see that, in general, 
the agreement is quite good, with the Eddington albedo being slightly 
lower than that given by the two-stream solution in all cases. For an 
excellent discussion of the relative accuracy of these two methods, com- 
pared with some exact results, see the review article by Irvine and 
Lenoble (1973). 

The Delta-Eddington Method 

Possibly the most challenging problem to be faced in radiative trans- 
fer theory is how to handle very asymmetric phase functions. Neither 
the simple two-stream approximation nor the Eddington solution can 
adequately cope with a phase function with a sharply scattered for- 
ward peak, which is the type usually encountered in aerosol and cloud 
studies. 

The delta-Eddington method was devised by Joseph, Wiscombe, 
and Weinman (1976) to allow the computationally simple Eddington 
method to be applied to sharply peaked forward-scattered phase func- 
tions. They approximate the phase function by a Dirac delta to account 
for a portion of the forward peak, and a two- term expansion for the rest 
of the phase function. The two-term expansion of P{cos9) follows from 
equation (2-31) as 


P{cos0) ^ 1 + 0 JiPi{cos 6 ) 

and with the definition of the asymmetry factor g of equation (5-57) 
this becomes 

P(cos 6 ) = \-\- igP\ (cos 6 ) 

The delta-Eddington phase function can now be written as 

P(cos0) — 2/6(1 - COS0) -h (1 - /)(! -{- 3gcos0) (6-23) 
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Two-Stream 

Eddington 

• Exact (doubling) method 
o Del ta-Eddington 



Figure 6-1. Comparison among four methods for a semi-infinite atmosphere. The 
dark circles are for the exact (doubling) method. Henyey- Greens tein phase 
function with g = 0.5. 
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where / is the fraction of the radiation scattered into the forward peak. 
The azimuthal averaging of the phase function follows as 

1 

= — / P(cos0) d(f> (6-24) 

27t Jo 

and using equation (1-33) with 

5(1 — COS0) = 27r5(/i — p^)6{(j> — <j/) 
we write equation (6-24) as 

F(/i, /) = 2/<5(/x - n') + (1 - /)(1 + (6-25) 

NoWj from the azimuthally averaged form of the RTE, equa- 
tion (2-41) 

/■ 1 

=1 y 2/«(/x - /)/(r,/) d/x' 

+ J j (1 - /)(1 + 3s/x/x')/(r, /) V 

/•I 

=u)fI{T,ix) + J (1 + Zgiin')I{r,n') d/i’ 


which can be rearranged to give 

"(/-i/Hr “ = I /_, + 3!W‘/)nr,K') V (6-26) 

But this is precisely the form of the RTE, equation (6-2), if we define 

r' = (l-w/)r (6-27) 


and 


OJ 


I — (jjf 


(6-28) 


Joseph et al. also show that for the delta-Eddington phase function to 
have the same asymmetry factor as the original phase function (the one 
we are trying to approximate), then 


9 = f + {i-f)g' 
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or 



9- f 
1 -/ 


(6-29) 


Finally, if the Henyey-Greenstein phase function is used, which does 
indeed produce reasonable results for many applications, they show 


that 




(6-30) 


and hence, equation (6-29) becomes 


9 

1 + g 


(6-31) 


Thus, the same equations derived earlier for the Eddington solution 
may be used if we replace u), r, and g with Q', t', and g*. 


Comparison of Two-Stream and Eddington Results 

The two-stream, Eddington, and delta-Eddington solutions are 
compared with the exact (doubling) method (Liou, 1980) in figures 6-2 
to 6-5, for the case of conservative scattering (w = 1) and one case of 
nonconservative scattering (d) = 0.8) for two optical depths, r* = 4.0 
and T* = 0.25, with g = 0.25 used in both cases. The superiority of 
the delta-Eddington method is clearly evident, especially at the more 
nearly vertical incident angles [p ^ 1.0). 
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Figure 6-2. Comparison of the two-stream, Eddington, and delta-Eddington methods for a finite atmosphere. The dark circles 
are for the exact (doubling) method. (Henyey-Greenstein phase function with r* = 4.0 and g = 0.25.) 
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Discrete Ordinates Method 

The discrete ordinates method is a very powerful analytic approach 
to solving the RTE. This method was developed by, or at least perfected 
and popularized by, Chandrasekhar. The procedure can be used to ex- 
tract numerical results for the simpler forms of the phase function, and 
has been used for numerical studies of nonhomogeneous atmospheres. 

Its greatest utility, however, seems to be a starting point for many 
theoretical attacks on the RTE. The theory has been developed to a very 
high degree of sophistication, and for that reason, it is worth spending 
some detailed effort in introducing this approach. The mathematics 
appears formidable at first glance, but once the reader gets into it, it 
emerges much simpler than imagined. 

The analysis here is confined to homogeneous semi-infinite atmo- 
spheres, and is carried far enough to permit the introduction of the 
well-known i7-function of Chandrasekhar. The principle of invariance 
for semi-infinite atmospheres is introduced and the integral-equation 
formulation of the iiT-function derived. The zeroth and first-order so- 
lutions to the integral equation are also derived, and some numerical 
results are given for higher-order approximations. Finally, some ele- 
mentary applications are presented. 

The extension to finite atmospheres is not given here, as it would 
be far beyond the intended scope of these notes. However, once the 
semi-infinite atmosphere case is understood, the reader will have little 
difficulty extrapolating to the finite atmosphere development, and the 
X- and F-functions, which are the finite atmosphere analogues to the 
H-function, will no longer seem quite so formidable or incomprehensi- 
ble. 

The analysis to be presented here follows chapter 3 of Chandrasekhar 
very closely, and merely supplies some of the missing steps in his 
development, although his text is so well written that it is difficult, 
even within the confines of the present supplement, to improve on it 
much. The analysis is restricted to two cases; (1) conservative isotropic 
scattering, and (2) nonconservative isotropic scattering. Again, his 
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results for nonisotropic scattering can easily be followed once the 
isotropic case is understood. 

To repeat something stated earlier, the isotropic case should not be 
dismissed lightly. As discussed in Irvine and Lenoble (1973) and Sobolev 
(1975), it is possible to develop similarity relations which can be used 
in some cases to approximately reduce an anisotropic scatter problem 
to an equivalent isotropic one. These relations allow an equivalent 
isotropic optical depth and single-scatter albedo to be defined in terms 
of the real anisotropic parameters. The isotropic problem is then 
solved and the solution transformed back to the “real” problem space. 
Similarity relations will be discussed briefly at the end of the chapter. 


Gaussian Integration 


We first present a few identities derived from Gaussian integration, 
as some of these results are needed in later developments. 

Basically, the integral of a continuous function is replaced by a finite 
sum 



m 

j~—m 


(7-1) 


where the weights Oj are given by 


1 


ay — 




n Pmjx) 
J — 1 X X y 


dx 


(7-2) 


and the ordinates xj are the zeros of the Legendre polynomials, P-mix). 
For our present needs, it is convenient to restrict ourselves to the zeros 
of the even-numbered polynomials, P2m(^)* (See the discussion in 
Chandrasekhar for more details as to why this is so.) For these divisions 


Q^i a — itX — i — 

For the case in which / (x) = x^ we get 

2 


(7-3) 


/: 


x^ dx = (for m even) 

m -h 1 ^ ^ 

— 0 (for m odd) 


Then, since 


p^x^dx= 


J=-n 
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E 

j = -n 


„ jm ^ 

m + 1 


= 0 


(for m even) 
(for m odd) 


(7-4) 


Abramowitz and Stegun (1970) give tables of aj and Xj for a number of 
orders n. See Chandrasekhar, or any good text on numerical methods, 
for more details of the Gaussian method. 


RTE for Conservative, Isotropic Scattering 

The governing equation for this problem is equation (2-41), with 
Cj — 1 and P(cos 6) — I 

= /(r, ^l)-\ f ^ I{r, mO V (7-5) 

Replace the integral with the Gaussian approximation and evaluate 
equation (7-5) at each of the 2n streams defined by the Gaussian 
quadrature points (see fig. 7-1). 


Mi 


dr 


1 ” 

= ^ ttjlj 


(7-6) 


Thus, equation (7-6) becomes a system of 2n linear equations with 
constant coefficients. As usual with systems of differential equations of 
this type, assume a set of exponential solutions 

I^=g^e-^^ (f=l,2,...,n) (7-7) 

where the are unknown constant coefficients. Substitute equa- 
tion (7-7) into equation (7-6) and reduce 

1 ^ 

ffi(H-Mifc) = 2 ^ 

j=-n 


Now, even though we do not know what the numerical values of the 
Qi are, they are constants, and the right-hand side of equation (7-8) 
directs us to sum over all these constants. Thus, the right-hand side of 
equation (7-8) is also a constant, K, and therefore, the must be of 


the form 


K 

~ 1 + mk 


(7-9) 
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Figure 7-1. Sketch showing the ray directions for the n = 4 case discussed in the 
text as the running numerical example. 


If we substitute equation (7-9) back into equation (7-8), we get 


1 

2 . 


E 


J = — 71 


Oj 


1 + Hjk 


= 1 


(7-10) 


Now, the limits on the sum are from —n to We can use equa- 
tion (7-3) to simplify equation (7-10) and write it in a neater form. If 
we expand equation (7-10) and look at the j = m term 


1 1 j 

l + l-pimk J 
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Since each j produces a pair of terms like this, we can write, using 
equation (7-3) 

2ctm 


1 

^ = 2 


+ 


1 - 

SO that we can write equation (7-10) as 


This is the characteristic equation for the equation set (7-6), from which 
we can get the 2n eigenvalues, Equation (7-11) is of degree n in 
and, thus, it can be seen that the eigenvalues occur in pairs, ±fca- For 
= 0, we have from equation (7-11) 

n 

E 1 

i=i 


while equation (7-4) gives, for m = 0 

n n 

E ai = 2^^ay = l 

j=-n j=l 

and hence, == 0, or A: — 0 is also a double root of equation (7-11). 
Note that this results from the assumed conservative scattering. 

Note that equation (7-11) has n vertical asymptotes (see, for ex- 
ample, fig. 7-2) — namely those values that occur at fc = 1/Pi> If we 
write 

^ ay 

we see that F(0) = 0, and for fc 1, F{lc) > 0. Also, 

lim F(k) = -1 
fc — 0 

Thus, F{k) plots as shown in figure 7-2, where we have used n = 4 as 
an example. (The 4-point Gaussian example will be carried throughout 
this chapter.) Since /i(= cos 9) < 1, the eigenvalues are positive, with 
one root at A: = 0. The roots can be found by the Newton-Raphson 
method ipn. 
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where the starting values can be taken as {l/p + c), where e is some 
small number. 

For n = 4 (4-point Gaussian quadrature), Abramowitz and Stegun 
(1970) give the following; 


= 0.1834346425 
P2 = 0.5255324099 
(is = 0.7966664774 
^4 = 0.9602898565 


ai = 0.3626837834 

02 = 0.3137066459 

03 = 0.2223810345 

04 = 0.1012285363 


Using equation (7-12), we get the roots in table 7-1. 


TABLE 7-1. ROOTS FROM EQUATION (7-12) 


a 

ka 

0 

0. 

1 

1.103185321 

2 

1.591778876 

3 

4.458085719 


With the ka now given, equation (7-7) gives the complete solution 
for as 


71 — 1 


a=l 


y 1 + j 


n— 1 

+ E 


a=l 


V 1 - t^i^a j 


(7-13) 


where and L-g are constants of integration. 

But, we have not yet included the k = 0 root. Guided by the grey 
Eddington solution (see e.g., Kourganoff, 1963), we assume for I( the 
solution 

I- = 6(r + qi) (7-14) 

with 6 and constants. Substitute equation (7-14) into equation (7-6) 



n 


j--n 


(7-15) 


and these equations can be satisfied if we let 

Qi — Q Pi 


(7-16) 
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where Q is another constant. Thus, with equation (7-13), the complete 
solution to equation (7-5) is 



+ T + Pi + Q 


(7-17) 


In equation (7-17), the b, Q, and L±ct {oc = 1, 2, . . . , n — 1) are the 2n 
constants of integration. 

We can eliminate some of these immediately. The radiance should 
not, of course, become infinite as r ^ oo. Thus, all the L_q. must 
vanish, and equation (7-17) reduces to 


h 


= b 


~n—l 

E 


La=l 


1 -I- Pika ) 


+ T + Pi + Q 


(7-18) 


One relation among the remaining constants can be found by applying 
the boundary condition that the incoming diffuse radiation at the top 
of the atmosphere (r = 0) be zero for all the -ju,-. This gives, for 
equation (7-18), 


n— 1 


o=Ei 


Qf=l 


^ikc 


-fli -\- Q 


(7-19) 


Equation (7-19) gives n equations in n unknowns, Q and n - 1 values 
of La> The constant b is not as yet found — it is left arbitrary for now. 
Thus, of the 2n original constants of integration, n — 1 are found to be 
zero by the requirement that remain finite as r oc, n are found 
from equation (7-19), and 6 is as yet unknown. 

Somewhat later in his text, Chandrasekhar goes to great lengths 
to develop a direct and simple way to determine numerical values for 
the constants La and Q — probably because at the time his original 
text was written, there were no efficient and accurate methods for 
inverting large matrices or for solving large systems of linear equations. 
With modern computers and numerical techniques, these sophisticated 
algebraic methods are no longer needed and will not be developed 
in these notes. Instead, we will solve directly the system of equa- 
tion (7-19). 
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For the 4-point Gaussian example, equation (7-19) gives the system 
of equations 

1.253703Li -I- 1.412414L2 -h 5.487454^3 + Q = 0.1834346 
2.379598Li -h 6.117365L2 - 0.744676L3 + Q = 0.5255324 
8.255792Li - 3.729726^2 - 0.391910L3 -|- Q = 0.7966655 
-16.840604Li - 1.891903L2 - 0.304781L3 + Q = 0.9602899 

The solution to these equations gives 

Li = -0.009461126 
L2 -= -0.036186730 
Ls = -0.083921097 
Q= 0.706919484 

(The values in Chandrasekhar are inadvertently given in reverse order. 
A note found in Kourganoff, 1963, p. 104, points out this reversal.) 


Some Elementary Identities 

Note that the solution of equation (7-18) contains a term similar to 


n— 1 

E 


Q=1 


1 

1 -h fiika 


It will be convenient to generalize this to a continuous function, and 
define the moments 

t 

We can derive a recursion formula for Dm{x) 


D 


m— 1 


w = E 


1 -h ^iX 


=E“.'‘r‘(' 








_fH^\ 

1+PixJ 
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But by equation (7-4) 




771— 1 

i 


hence 


or 


Dr 


m 


2(5 


_ 1 W = - - X - rtD„(x) 

m “ 1 + u,x m 


(<5 = 0 for m even) 
(^ = 1 for m odd) 

26 


'^^777(^) — ~ 
X 


-h p^x m 


2^ n f ^ 

1 (^) 
m 


is the required recursion equation. Thus we get, for m odd 




1 


x [2j 


- - D2j-2{x) 


and for m even 


-1 




(7-21) 

(7-22) 

(7-23) 


By comparing equation (7-20), with m = 0, to equation (7-10) 

Oo(-,)= = 2 (^24) 

and, thus, from equations (7-22) and (7-23) we get the sequence 

Do (a;) = 2 
Di(x) = 0 

D2{x) = 0 

= i 


By repeated application of equations (7-22) and (7-23), with the above, 
it is possible to establish the general formulas (see Chandrasekhar), 
since A: is a root of the characteristic equation. 


D2j-l{k) = 


+ 


[2j -l)k {2j - 3)fc2 




3A;2a-3 


( 7 - 25 ) 
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-2 


" (2j - l)fc2 (2i - 3)fc4 ■ ■ ■ 3fc2j-2 


(7-26) 


wherej = 2,3,...,n. . r 

The even Legendre polynomials can always be written m the torm 


P2nilA = 

J=0 


where the P 2 j are constant coefficients. 
Consider the expansion 


ajpf 


E P2jD2j{k) = E P2i E r+ mfc 

J=0 ^=0 




j=-n 
n 


i—0 

m 


- E rr1^E/-w 


J--71 


z=0 


(7-27) 


But, by the Gaussian quadrature procedure we have adopted, the pj 
are zeros of the even-numbered Legendre polynomials. Hence, the right- 
hand side of the above equation is equal to zero, and we have 


Ep2,%W=0 (7-28) 

j=0 

From this Chandrasekhar derives an equation which is used in several 
places in the remaining development 

(fci, fc2 . • ■ , kn-i){pi, ■ 1 Mn) = (7-29) 

(Note the different ranges on the subscripts of equation (7-29).) 


The Flux Equation 

The flux is defined as before, for azimuthally symmetric radiation 
F[t) = p,I{T,fi)dp. 
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Or, in the Gaussian approximation 


n 

F(t) = 27t (7-30) 

j=-n 

If we substitute the solution (7-18) into equation (7-30) 


F{t) - 2irb 


n~l 


z 

,a~l 



1 + fiika 


n n 

+ (’• + (?) ^ aifn+ ^ ai^t^ 
i—~n i=~n 


From the definitions of Dm{^) in equation (7-20) this becomes 


F{t) = 27t6 


^ Ti n 

Z.aC Di(k(x) + + Q) 


,Qt= 1 


But Di{ka) ~ 0, and from equation (7-4) 






2 

3 


and thus we get 

(7-31) 

Since 6 is a constant, this equation says that the flux is constant at 
all r— which is indeed true for this problem. (Since we have considered 
the equilibrium problem of conservative scattering in a semi-infinite 
atmosphere, the net flux m must equal the net flux out, and the net flux 
is conserved at all altitudes because there is no absorption or emission.) 
Equation (7-31) allows us now to evaluate the constant b 

b= -ttF 
4 

and, hence, this establishes our final constant of integration. The 
solution of equation (7-18) becomes 
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1 4 - Pika J 

in which all the constants are now known. 


+ T + Hi + Q\ 


(7-32) 


The Source Function 

From differential equation (7-5), the source function for this problem 
is written in the Gaussian quadrature form 



£ 


I[t,h) dfi 




1 

2 X) 

n 


(7-33) 


Insert the solution equation (7-32) into equation (7-33) and proceed 
as in the F-integral; i.e., interchange orders of summation and use the 
Dm{x) definitions, and we find that J reduces to 


J = -ttF 
4 




(7-34) 


Following Chandrasekhar, define 

n— 1 

?(t) = ^ -I- Q 

a=i 

and we can write the source function in the Eddington form 

J = ^7rF[r -I- g(r)] 

Inserting our numerical values in equation (7-35), we get 


(7-35) 


(7-36) 


g(r) =0.706919 - 0.009461 exp(-1.103188r) - 0.036187 exp(-1.591778T) 
- 0.83921 exp(-4.45808r) 


See table 7-2. 

Given the source function, we can now use equations (3-5) and (3-6) 
to get the intensity at any t, h (assuming, of course, the same boundary 
conditions and symmetry in = — /x). Thus, we write equations (3-5) 
and (3-6) as 


m 
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I{r.-n)= r 

Jo /i 

Substitution of equation (7’34) into the above and integrating poses no 
major problems. The results are 


3 

^(^1 m) = I X/ 1 , .. + ^ + M + Q 


\a=l 


1 + kafi 


(7-37) 


= -ttF 
4 


E + - + ( c ? - /^)(1 - 

,a = l 


(7-38) 


TABLE 7-2. VALUES OF q{T) FROM 
EQUATION (7-35) 


T 

Q{t) 

0. 

0.577350 

0.1 

0.613849 

1.0 

0.695441 

3.0 

0.706268 

5.0 

0.706868 

10.0 

0.706919 

oo 

Q 


These are the final forms for the upward and downward radiance 
components in the nth approximation. Note that this is the intensity in 
direction (i based on a 2n-stream approximation for the source function. 

The Law of Darkening 

By putting r = 0 in equation (7-37), we get the angular distribution 
of the radiation emerging from the top of the atmosphere; i.e., the law 
of darkening or the limb darkening equation 

= + + (7-39) 

From our numerical example, we get 

_ 3 r -0.009461 0.036187 0.083921 , I 

7tF ~ 4 L 1 + 1.103188m 1 + 1 159178m ~ 1 + 4.458080m ^ ^ ^ 0.706919J 
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TABLE 7-3. VALUES FROM EQUATION (7-39) 




imF) 

/(0,l) 

0.0 

0.433013 

0.345082 

0.1 

0.531852 

0.423850 

0.2 

0.620516 

0.494509 

0.3 

0.704562 

0.561488 

0.4 

0.786070 

0.626444 

0.5 

0.866012 

0.690152 

0.6 

0.944910 

0.753029 

0.7 

1.023094 

0.815321 

0.8 

1.100699 

0.877182 

0.9 

1.177915 

0.938718 

1.0 

1.254812 

1.0 


See table 7-3. 

Compare equation (7-18), a set of equations with a discrete argu- 
ment, Hi, with the parenthetically enclosed term of equation (7-39), 
which is a continuous function of p. These are identical in form, except 
for the sign of the //. As a lead-in to the ff-functions, we define the 
continuous function 


n— 1 

5(/i) - E 




1 — kafJ' 


— p-\- Q 


(7-40) 


From this, we can write the boundary conditions (7-19) as 

S{ni) = 0 (f = l,2,...,n) (7-41) 

and the law of darkening, equation (7-39), can be put into the form 




(7-42) 


a form which will be found to be more useful after we have derived the 
H-function. Note that the quantity 37r5(-/i)/4 can be to sonae extent 
interpreted from equation (7-42) as a diffuse reflection coefficient, i.e., 
it gives the angular distribution of the reflected radiance in terms of 
the incoming flux, F. 
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The //-Functions 

Note from the definition of equation (7-40) that the 5-function is 
defined in terms of the La and Q. These must in turn be obtained 
by solving a set of linear equations, equation (7-19), Knowledge of the 
La and Q permits us to determine the intensity and fiux components 
at any point within the atmosphere. However, in many cases we 
are not concerned about the detailed structure of the radiation field 
inside the medium, but really need to know only what comes out 
of the top and/or what comes out of the bottom of the atmosphere. 
Chandrasekhar presents a method for doing just this, and this analysis 
leads to the definition of the //-functions, a set of functions which, for a 
given phase function, can be computed once and for all and tabulated. 
Note carefully the distinction between emission and scattering in the 
following development. Equation (7-51) expresses emission in terms of 
the //-function, while equation (7-84) expresses scattering in terms of 
the //-function with two different arguments. We proceed now with 
this derivation. 

The summation in equation (7-40) contains the expression {l~kapf 
in the denominator of each term. If we define the function 


n— 1 

^(m) = rt ~ (7-43) 

a=l 

then, by multiplying S[p) by /?(//), we get a function which is clear of 
fractions. 


n~\ 

s{^l)R{n ) = n (^ ~ 

a=l 



L 


a 


1 — kaIJi 



(7-44) 


Since R{p) is a polynomial of degree (n — 1) in //, the presence of the 
//-term in the parentheses of equation (7-44) means that the product 
5{//) /?(//) is a polynomial of degree n in //. Also notice that 5(//)/?(/i) 
vanishes for // = //^, / = 1, 2, . . . , n, since by equation (7-41) 5 vanishes 
for these values. 

Define the polynomial 


i=\ 


(7-45) 


120 


Chapter 7 


which is also a polynomial of degree n in p. Since P{p) and S{p)R{p) 
are both of degree n in p, and have the same roots, pi, they can differ 
from each other by, at most, a multiplicative constant; i.e., 

S{p)R{p) = KP{p) (7-46) 

The constant K can be determined by comparing any power of p on 
both sides of equation (7-46). In this case, it is easiest to compare 
coefficients of the highest power; i.e., of p”'. From equation (7-44), it 
can be seen that the coefficient of p^ on the left-hand side is 

(-l)”/cifc2---fcn-l 


while the coefficient of p'^ on the right-hand side is unity. Thus 
and hence, from equation (7-46) 


S{p) = {-irkik2...kn-i^ (7-47) 

With the definition equations (7-43) and (7-45), and equation (7-29), 
this can all be put into the form (changing the sign of p) 


S{-p) 


1 1 
s/Z PlP2 ■ ■ ■ tJ-n 




1=1 
n— 1 


(1 -I- fca/i) 


a=l 


(7-48) 


which now only contains the discrete coordinates, /Zj, and the eigen- 
values, fca, of the original system of equations. From equation (7-48) 
comes the discrete form of the definition of the i^-function 




1 


+ Mi) 




MlM2 ---Mn 


n=l 

(1 -I- kaP) 


a=l 


(7-49) 
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The equation (7-48) can be written 


and the law of darkening, equation (7-42), becomes 


/(0,/i) = ^7TFi/(/z) 


(7-50) 


(7-51) 


For our numerical example, in which n = 4, we get from equation (7-49) 

1 


//(//) = 

^ 0.073749 

See table 7-4. 


(M + 0.18343)(// + 0.52553) + 0. 79667) (/i + 0.96029) 1 

(1 + 1.10319/i)(l + 1.59178//)(1 + 4.45808/i) J 


TABLE 7-4. VALUES FROM EQUATION (7-49) 



H {p)fi=4 

^(/^) exact 

0.0 

1.000000 

1.0000 

0.1 

1.228240 

1.2474 

0.2 

1.433003 

1.4503 

0.3 

1.627100 

1.6425 

0.4 

1.815335 

1.8293 

0.5 

1.999953 

2.0128 

0.6 

2.182162 

2.1941 

0.7 

2.362674 

2.3740 

0.8 

2.541940 

2.5527 

0.9 

2.720262 

2.7306 

1.0 

2.897849 

2.9078 


The exact solutions were computed from an integral equation to be 
developed later. It can be seen, however, that the fourth-order solution 
is not too bad, considering the relatively simple arithmetic involved. 

As mentioned earlier, Chandrasekhar goes to great lengths to de- 
velop expressions similar to equation (7-47) for the constants La and 
Q. These derivations will not be repeated here, for reasons already 
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stated, but his results will be given for completeness and reference. 
Define the polynomial 




n—1 


n (1 - 

0=1 

0^a 


(7-52) 


Then the constants La and Q can be found without having to invert 
any matrices, from the relations 


La 


. fcn-l 


P{l/kg) 
Pa ( 1 / ka ) 


(7-53) 


with P{x) given by equation (7-45), and 


n 71—1 ^ 

1 ^ 1 

i=i a=l 


(7-54) 


From our numerical example, we get 


P(l/fci) = -0.00162775 
P(l/fc2) = 0.00255488 
P(l/fc3) = -0.00518676 
Pi(l/A:i) = 1.346865 
P2(1/A;2) - -0.552715 
P3(1/A:3) =0.483843 

We also get kik 2 k^ - 7.828524. Then, using equations (7-53) and 
(7-54) we get 

Li = -0.0094611411 
L2 = -0.0361860937 
Ls = -0.0839211267 
Q= 0.70691923070 

which can be compared with the values found earlier, following equa- 
tion (7-19), by inverting a 4- by 4-matrix. 

Chandrasekhar also derives an accuracy check 


n—1 

g+ 


a=l 


J_ 


(7-55) 
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and, using our numerical data, we get for the left-hand side 

0. 5773508692. compared with the exact value of 0.5773502692. 

So far, we have achieved a number of significant goals: 

1. We have completely solved the simple case of the radiant field 
for which the net flux is constant, and a conservative, isotropic 
scattering medium. 

2. We have gotten some numerical values — in the fourth-order 
approximation — for some of our expressions. These do not seem 
quite so frightening any more. 

3. We have developed some basic concepts and ideas which will be of 
more use later on— specifically the S- and i7-functions. 

4. We have solved the limb darkening case for this simple problem. 

Now, we consider a somewhat more difficult and useful problem. 

Diffuse Radiation With Non-Conservative, Isotropic Scattering 

We will now consider the problem of the scattering of a collimated 
beam into a semi-infinite atmosphere; i,e., the scattering of sunlight 
by a planetary atmosphere. Our starting equation is equation (2-50), 
which, for isotropic scattering, becomes 

- I ( 7 _ 56 ) 

Here, u is the single-scattering albedo, and it is assumed that a parallel 
beam of solar radiation of flux F is hitting the top of the atmosphere 
at the angle Oq — cos~^ /iO- 

The development given here will be somewhat sketchier than that 
in the earlier part of this chapter, as they are very similar and should 
not now pose any difficulty. 

Again we discretize the integral and solve equation (7-56) along the 
discrete rays defined by 

n 

(i = ±1. ±2, . . ,±n) (7-57) 

j=-n 

We first solve the homogeneous system 

dl 1 ^ 

- 2^^ E (7-58) 

j=-n 


124 


Chapter 7 


by assuming, as before, 


If we put this into equation (7-58) and reduce, we get the characteristic 
equation for this problem 

This is identical to equation (7-11), except for the presence of the Cj, 
But this is a big “except,” for now there are no roots at fc = 0, and it 
will not now be necessary to introduce the somewhat artificial solution 
(7-14) into the system. The complete solution follows directly. 

The asymptotes of equation (7-59) occur at the same place as those 
of equation (7-11), but the roots are somewhat larger, depending on 
the value of (D, since 




i=i 




OJ 


> 1 


The roots of equation (7-59) can be found for any oj in the same way 
as before. For the n = 4 case the data are as in table 7-5. 


TABLE 7-5. ROOTS OF EQUATION (7-59) 


d) 

a — 1 

a = 2 

a = 3 

a = 4 

1.0 

0. 

1.103186 

1.591779 

4.458086 

0.9 

0.525430 

1.108937 

1.615640 

4.554851 

0.8 

0.710413 

1.116799 

1.642629 

4.652965 

0.7 

0.828671 

1.127655 

1.672473 

4.752078 

0.6 

0.907693 

1.142395 

1.704602 

4.851871 

0.5 

0.959481 

1.160900 

1.738275 

4.952060 

0.4 

0.992327 

1.181880 

1.772515 

5.052401 

0.3 

1.012963 

1.203057 

1.806656 

5.152683 

0.2 

1.026230 

1.222732 

1.840027 

5.252727 

0.1 

1.035120 

1.240155 

1.872179 

5.352384 
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So, a set of solutions to the homogeneous equations is 








(7-60) 


Now, we need a particular solution to equation (7-57). Assume one 
of the form 

h = ^QFhie-^/f^o (7_61) 

where the h^ are constants. Substitute into equation (7-57) and we find 
that the hi must satisfy 

'‘■(‘ + w) = 5" £ “A + ' 

j=~n 

and hence the hi must have the form 


“ 1 -t- ^ 
MO 


(7-63) 


where 7 is an unknown constant. Put this back into equation (7-62) 


7 


• (jj 


E 

j=i 




(7-64) 


Put equation (7-63) into equation (7-61) and combine with equa- 
tion (7-60), and we get the complete solution to the system of equa- 
tions (7-57) 




LL 


o-kaT 


L' 






76 


-'''/mo 


a=l 


- H — (jjF- 
fl^ka 4 1-1- ifli/flo) 


(7-65) 


As before, in order to bound the radiance as r ^ oo, we must require 
that all the L'_^ = 0, which leaves 






LctS 'ye ^0 

\+ ^lika 1 + {Fi/Fo) 


(7-66) 
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We apply the same boundary conditions at the top of the atmosphere — 
the incoming diffuse radiation at r = 0 is zero along the rays — and 
we get the system of equations for the n remaining constants La- 




+ 


a=l 


l^i^a 1 (Mi/mo) 


= 0 


(7-67) 


As a matter of comparison, note the difference between equation (7-67) 
and equation (7-19) for the conservative scattering case. In equa- 
tion (7-67) the a-summation goes from 1 to n rather than from 1 
to (n - 1), and there is no Q constant. As pointed out above, these 
differences result from the fact that = 0 is not a root. All the 
solutions we need are contained in equation (7-67). 

With 7 defined by equation (7-64), equation (7-67) again provides 
us with n equations in the n unknowns La- We can solve this system 
just as before to get the complete solution for the total radiation field. 
However, if we only want the law of darkening for the emerging field at 
the top of the atmosphere, this can again be expressed in terms of the 
77-function, and we need not evaluate the La- However, to carry along 
the numerical example, we will evaluate equation (7-67) for the n = 4 
case we have been using. We will put arbitrarily u> = 0.8 and get 

7 = 0.947722 
Li =0.516131 
L2 = 0.046078 
L3 = 0.246943 
L4 = -0.402956 

This source function for this problem in the Gaussian approximation 
is, from equations (7-57), 

Ar) = li E «>/, + 

j=-n 

If we substitute equation (7-65) into the above and reduce, 

J(r) = ^u)F LaC~'^°‘^ -h 76“’’/'^°^ 
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If we substitute this into equations (3-5) and (3-6) as before, with zero 
boundary conditions on the incoming diffuse radiation, we get for the 
complete solution to the radiation field at any r, fi 


(7-68) 


4 


^ l-kau'- ' I- u/m ' ’ 

,a = l 


The law of darkening follows from equation (7-68) with r = 0 


(7-69) 




^a=l 


-I- kaP P0 + 


(7-70) 


We would like to write equation (7-70) in terms of the i7-function as 
we did before. Again guided by the form of the characteristic equation, 
putting 2 = 1/fc we write for equation (7-59) 


1 


1 - 



and define the continuous function 


T(^) = l-(Z;^2^ 
J=1 



(7-71) 


(7-72) 


Obviously, this must vanish for ^ = l/k, since A: is a root of the 
characteristic equation. 

Now, T{z) is a polynomial of degree zero. Thus, 


nz)i[{z‘^~pj) 

is a polynomial of degree n in 2 with roots ±l/ka, a = 1,2,. . . ,n. The 
polynomial 

J = 1 
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is also of degree n in z with roots ifcl/fca- Thus, as before, these two 
polynomials can differ by at most a multiplicative constant 

fi ^ n (i - 

j=i j=i 

From equation (7-72) we see that T(0) = 1, and hence if we set z = 0 
in equation (7-73) 

K = (-l)Vf/x| -/i2 

and, thus, from equation (7-73) we can write T{z) as 






01 = 1 


(2^ - /*«) 

n n 

(1 - kaz) 1 ^ (1 + 




(7-74) 


(- 1)" n n 


Of=l Q=1 


But from the defining equation (7-49), this can be written 

^ H{z)H{-z) 

If we let z = po in equation (7-72) 


(7-75) 




T(/^o) — 1 ^ 2 ,,2 

j=l ^0 ^ 




y=i 


But this is exactly equal to the denominator of equation (7-64); thus, 
we get 

7 = (7-76) 

giving us the unknown constant 7 in terms of the 77-functions. 
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Now, again repeating the earlier procedure, we can be guided by the 
form of the law of darkening, equation (7-70), and with equation (7-76) 
define the continuous function 




a— 1 


1 - (m/mo) 


(7-77) 


in which again S{pi) = 0, f = 1, 2, . . . , n. The law of darkening becomes 


7T , 




(7-78) 


which we want to write in terms of the i7-functions. The function 

n 


(l-f)5(//)Il(l-W) 

^0 ail 


(7-79) 


is a polynomial of degree n in //, which vanishes for (x = pi because 
S{Pi) = 0. Thus we can write 


Hl 

MO 


5( m ) (7-80) 

a=l a=l 


or 


, n 

a=l 

Comparison of equation (7-81) with equation (7-48) shows that this is 
almost in the right form. If we redefine the constant K' 


K' = K 


(-ly 


• • Mn 

then equation (7-81) can be manipulated to the form 

1 - 

in which we now need to evaluate K. From equation (7-77) 


(7-82) 


1 


-) 

ml 


5(m) 
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and from this we can see that as p approaches po 

lim (l- S{p) = H{po)H{-po) 
^^^0 V MO/ 

while from equation (7-82) 


lim (l--^]s{p) = KH{-po) 

\ pqJ 

and, thus, we find that K = H{po), and equation (7-82) becomes 


S[p) = 


H{po)H{-p) 

1 - 


and the law of darkening, equation (7-78), becomes 

7(0, m) = "^uF-^H[po)H{p) 

4 po + p 


(7-83) 


(7-84) 


Recall Chandrasekhar’s definition of the scattering function, equa- 
tion (4-4) 

7T 

7(0, m) = 

(The 7 T comes in because of the way we have defined the incoming solar 
fiux— Chandrasekhar defines it to be ttF, while we have defined it as 
just F.) 

If we compare equation (7-85) with equation (7-84), we find 

l^F-f^H{po)H{p) = j-S{p,po) 

4 P0 + p 4 p 


or 

(- + —] S{p, po) = CjH{po)H{p) (7-86) 

Vm mo/ 

which gives the scattering function in terms of the tabulated 77- 
functions. 

Use of equation (4-10) allows the reflection function F(/i,Mo) also 
to be written in terms of the //-function 

(i + -) «,) = (7-87) 

\p PoJ 4/x/io 
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We note that interchanging p and pQ in equations (7-86) and (7-87) 
gives 

= S{pa,p) 

^(m,Mo) = R{(iQ,ij) 

which are examples of the law of reciprocity, a concept which occurs 
frequently and is much used in theoretical analyses. 

Applying our numerical example to equation (7-70), with po = 0.4, 
u) = 0.8, and n = 4, we get the following comparison in table 7-6 with 
Chandrasekhar’s exact results for the reflection function. Again we note 
that the n = 4 approximation is reasonably good, the maximum error 
being about 1 percent aX p — 0.5. 

TABLE 7-6. VALUES FROM EQUATION (7-70) 



II 

V /exax:t 

0.0 

0.270783 

0.272217 

,1 

.243725 

.248014 

.2 

.219711 

.222972 

.3 

.199753 

.202310 

.4 

.183164 

.185255 

.5 

.169225 

.170984 

.6 

.157331 

.158864 

.7 

.147081 

.148436 

.8 

.138144 

.139359 

.9 

.130277 

.131380 

1.0 

.123693 

.124303 


Similarity Relations 

We can note from our earlier developments that the radiation field is 
essentially characterized by three basic quantities — the phase function, 
the single-scattering albedo, and the optical thickness. A change in 
any of these quantities produces a change in the radiation field. The 
question thus arises: Can we change these parameters simultaneously 
in such a way that the radiation field remains at least approximately 
fixed? In particular, can we relate a given set of parameters P{p, p'), Co, 
and r to an equivalent set of isotropic parameters? 

The answer is obviously yes, or we would not have raised the question 
here, and this section would not have been written. We have already 
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seen one example of such transformations in the discussion of the delta- 
Eddington method of chapter 6, in that equations (6-27) and (6-28) 
give a transformation which relates the delta-Eddington solution to the 
classical Eddington solution. 

As pointed out by Sobolev (1975), the approximate similarity of 
the radiation field in an atmosphere with anisotropic scattering to 
the corresponding field in an atmosphere with isotropic scattering 
will take place only after a large number of scatterings; i.e., large 
optical thickness and u ^ I, Also, similarity can only be discussed in 
connection with azimuthally averaged fields, since the isotropic radiance 
is azimuthally independent. 

The diffuse radiation field in a plane-parallel atmosphere follows the 
now-familiar equation (2-29) 

_ il / I{t)P{cos0) dn (7-88) 

dr 47T J 

Suppose now we assume that the fraction r of the radiance is scattered 
isotropically (i.e., P(cos0) - 1), and the remainder (1 - r) is approx- 
imated by a Dirac delta function, so that we can write for the phase 
function 

P(cos 0) — r + (1 — r)6 (7-89) 

If we put equation (7-89) into equation (7-88) 


/xf = /-|^//(r)[r-h(l-r)^]dn 

= I-—r[ I{t) dU - ^(1 - r)47T7 
47T j 47T 

= /[/ - w(l - r)] j I{t) dn 


If we divide through by [/ — di(l — r)] 

a — ^ = / - — r / I{t) dn 

^ [/ — w(l - r)] dr 47T 1 — w(l — r) J 

and thus, by defining 


r/ = [1 -w(l -r)]r 


and 


1— o;(l — r) 


(7-90) 


(7-91) 


(7-92) 

(7-93) 


133 



Introduction to the Theory of Atmospheric Radiative Transfer 
we can write equation (7-91) in the form 



(7-94) 


Equation (7-94) is therefore identical in form to equation (7-88) with 
P{cos6) = 1, which thus describes isotropic scattering. In this way, 
equations (7-92) and (7-93) can be considered to be similarity relations 
which transform the anisotropic problem of equation (7-88) to the 
equivalent isotropic problem of equation (7-94), under the assumption 
in equation (7-89). 

We now have to determine the quantity r in equation (7-89). The 
more forward scattering we have, the smaller the value of r. We have 
seen earlier that, in the Henyey-Greenstein phase function, the factor g 
controls the size of the forward-scattering peak; the larger the amount 
of forward scatter the more nearly g approached unity. Thus, if we 
choose 

r - 1 - g (7-95) 

then we get 

ri - {l-ug)r (7-96) 


and 


OJl 


- g) 
I -u>g 


(7-97) 


as our set of similarity relations. 

The discussions in Sobolev (1975) and Irvine (1975) indicate that 
equations (7-96) and (7-97) produce solutions that agree well with more 
nearly exact solutions in most cases, the agreement generally being 
better for integrated quantities, such as total albedo or the atmospheric 
flux, rather than quantities such as radiance. Again, this is because the 
integration, after multiple scattering, tends to smooth out the effects of 
the phase function. 

The similarity solutions for one of the cases given earlier for the two- 
stream and Eddington methods are presented below and in figure 7-3, 
The were interpolated in the tables of /f-functions at the end 

of chapter 8. The similarity relations give the correspondences in 
table 7-7. The correspondences in table 7-7 in turn give the values 
for the reflection function shown in table 7-8. 
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TABLE 7-7. CORRESPONDENCES 


OJ 

(jJ2 

0.99 

0.9802 

.95 

.9048 

.90 

.8182 


TABLE 7-8. r(/xo), SIMILARITY SOLUTIONS 


MO 

u = 0.99 

Cj = 0.95 

d ) — 0.90 

0.0 

0.8593 

0.6915 

0.5736 

,1 

.8285 

.6377 

.5121 

.2 

.8048 

.6003 

.4718 

.3 

.7835 

.5688 

.4392 

.4 

.7635 

.5414 

.4116 

.5 

.7447 

.5169 

.3877 

.6 

.7268 

.4947 

.3667 

.7 

.7098 

.4746 

.3479 

.8 

.6971 

.4561 

.3311 

.9 

.6772 

.4391 

.3159 

1.0 

.6616 

.4234 

.3020 
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The Principle of Invariance 

The principle of invariance is a very elegant concept that was 
first put forth by Ambartsumyan (1958) and later perfected by Chan- 
drasekhar and others. The principle is deceptively simple and will be 
applied here to the case of a semi-infinite, homogeneous atmosphere. It 
allows us to derive a single integral equation for a function which can be 
identified with the i/-function. This linear equation permits the exact 
computation of the i?-function. (The numerical solution is an iterative 
one. Thus, it is exact only in the limit, but practically converges to 6 
to 8 decimals in a few iterations for small values of the single-scattering 
albedo. The convergence is slower as the single-scattering albedo ap>- 
proaches unity.) 

The principle of invariance for an infinitely thick atmosphere can 
be stated as follows: we are given the infinitely thick atmosphere with 
certain reflection and absorption properties. If we add an additional 
layer of the same optical properties to the top of the atmosphere, we do 
not change the overall reflective and absorptance characteristics of the 
atmosphere. By adding a thin layer, we can compute the differential 
change in reflection and absorption and set these changes to zero. The 
result is a linear integral equation for a function which we can relate to 
the /f-function derived in the last chapter. 

We will follow essentially the development of Liou, and use his 
definition of the reflection and transmission functions, and then relate 
the final equation to the form developed by Chandrasekhar. 

We assume that the added layer is so thin that at most a single 
scatter can occur in it. Then, for a given photon which is reflected out 
of the top of the atmosphere, only one of the five histories sketched in 
figure 8-1 can occur: 

1. The photon can penetrate the thin layer and be reflected from the 
infinitely thick layer (ITL). 

2. The photon can be singly scattered upward from the thin layer 
before it reaches the ITL, 
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Figure 8-1. Sketch showing the five single- scat ter scenarios between the added 
thin layer and the infinitely thick layer (ITL). 

3. The photon can be singly scattered downward by the thin layer, 
then reflect upward from the ITL. 

4. The photon can penetrate the thin layer, reflect from the ITL, and 
then be singly scattered upward by the thin layer. 

5. The photon can penetrate the thin layer, reflect upward from the 
ITL, and then be singly scattered downward by the thin layer to be 
once more reflected up and out by the ITL. 

We assume that Ar, the thickness of the thin layer, is 1, and hence 
only terms linear in Ar will be retained. The single-scattering albedo 
determines the fraction of the incoming photons which are scattered. 
In the thin layer, it is assumed that there is no absorption along a path 
which involves a single scatter, but that there is absorption along a path 
which penetrates the thin layer and along which there is no scatter. In 
other words, we assume that a given photon may be either absorbed or 
scattered, but not both. 

For azimuthal symmetry, Liou’s definition for the reflection function 
follows from equation (4-6) 

I{0, m) = 2 / Rifi, n')I{0, -n')n' dfi' (8-1) 

Jo 

The reflection coefficient for an infinitely thin layer can be obtained 
from the single-scattering solution given by equation (5-33b), which for 
r* = Ar 1 reduces to 

^( 0 > -fio)f^oFo ( 8 - 2 ) 

4fino 

and hence, the reflection coefficient becomes 

mo) == -//o) (8-3) 
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For thin layers, the transmission function reduces to 

g-Ar/M R, 1 _ 

The simplest way the writer has found to derive the differential 
changes in R due to the addition of the thin layer is to start with the 
emergent beam and work backwards to the source. This will be done in 
the five parts of figure 8-1 for each of the five scenarios sketched above. 

In figure 8- 1(a), reflection from the ITL, 

/(-Ar,^) = (1 - Ar//i)/(0,/x) 

/(0,/z) = 7(0, -no)R{fi,no) 

-/^o) = (1 - ^'r/no)noFo 



Figure 8- 1(a). Sketch of the first event, reflection from the ITL. 

Put all these together 

7(-Ar,/z) = (1 - ~ AT//io)jUo^b 

Expand and retain only terms to first order in Ar 

I{-At,h)/hqFq = R{n,fio) - R{n,fio){AT/fi + Ar/fio) 

But this is the new reflection coefficient, and hence the change in the 
reflection coefficient due to this first event is 

ARi{(i,fio) = -R{n,fio) At {l/fiQ + l//x) (8-4) 

In figure 8- 1(b), single upward scatter from the thin layer, 

I{-At,h) = R{fi,fio)fioFo 
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T = 0 


Figure 8- 1(b). Sketch of the second event, single upward scatter from the thin 
layer. 



Figure 8-l(c). Sketch of the third event, single scatter from the thin layer followed 
by a reflection from the ITL. 

and from equation (8-3) 

= ‘^^-^P{n,-no)fioFo 
4nfio 


and hence, 




ui At 
4fi/io 


P{n, - m ) 


(8-5) 


In figure 8- 1(c), single scatter from the thin layer followed by a 
reflection from the ITL, 

I{-At,h) = (1 - Ar//i)/(0,/i) 
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But since all possible /x' must be included, we must use equa- 
tion (8-1) 

/(-Ar,/i) = 2^ n' d/j.' 

= (l - ^) ^ R{n,n')P{-fi', -no) dn' 

and so to order Ar 

A p J 

AR3{n,no)^^^ R{n.n')P{-f^',-fj-o) dn' (8-6) 

2^0 Jo 

In figure 8- 1(d), reflection from the ITL followed by an upward 
scatter from the thin layer, 

i{-AT,n) = (^) P{n,n')mt^') 

7(0, n') = -Mo)^(o, -mo) 

•^(0, -Mo) = (1 - At/mo)mo^O 

Ji) 

T = -At 


T = 0 

Figure 8- 1(d). Sketch of the fourth event, reflection from the ITL followed by an 
upward scatter from the thin layer. 

Again, using equation (8-1) we get 

7(-Ar, m) = ^^^-^noFo ( 1 - — ) f R{n', -Mo)^(m, m') dn' 

2m \ mo / Jo 
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or 

j V (8-7) 

2 /^ JQ 

In figure 8- 1(e), reflection from the ITL, followed by a downward 
scatter from the thin layer, and a final reflection from the ITL, 


I{-AT,fj.) = (1 - Ar//i)/(0,/x) 



1(0, y") 1(0, y') 


Figure 8-1 (e). Sketch of the fifth event, reflection from the ITL, followed by a 
downward scatter from the thin layer, and a final reflection from the ITL. 

and hence, integrating over all 

/(0,/i) ^2 f h'R{h,h')I{0,-h') dfj,' 

Jo 


Now, 

iM = 

and we can write 7(0, //) as 

7(0, /i) = f R{fi,fj/)P{-fi.',n") dfi' 

2// 7o 

and to order At 

/(-Ar,/i) = ^^^7(0,/') f R{fi,fi')P{-(x',n") d^' 
2fi" Jo 
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But now, 

7(0, /i") = 7(0, -//o)72(/i", ~pq) = f 1 - — ) -Mo)mo^O 

\ A^O / 

But this must be integrated once more, this time over all jj!' , and we 
get 


/(-At,m) = 2(moFo) / 

Jo Jo 

or, regrouping the integrals 

AR^{n,(io) = uj At [ R[h,r!) dn' f R{fj.",-no)P{-n',n") dfi" 
Jo Jo 

( 8 - 8 ) 

Now, according to the principle of invariance 


A7?i + A7?2 + A7?3 + A 7 Z 4 -j- A7?5 — 0 


(8-9) 


and so, from equation (8-4) to equation (8-8) in equation (8-9), we 
divide out the At and factor out d;/4/z/i0) write 

f — + =7^ +2/X / V 

( Jo 

) f dfi' 

Jo 

i I P(/i, //') rf/i' 

J2M0 J V'j I 


+ 2/iO 


+ 2^ 


(8-10) 


which is the desired integral equation for /?(//, /^o)* Note that this 
equation is nonlinear. The only restrictions on equation (8-10) are that 
the atmosphere must be homogeneous, plane-parallel, and semi-infinite. 

Now, let us consider the case of isotropic scattering. Then all of the 
phase functions in equation (8-10) are equal to unity, and 



P(/x,Mo) 


4///X0 
H- 


\~2p f R(p,p) dp ^2po ( R{p\po)dp 
Jo Jo 

4ppo J R{p,p)dpj P(/x",/[io) d/i"| (8-11) 
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If we interchange p and po equation (8-1 1), we get the same expres- 
sion, indicating that R{p,po) and R{po,p) both satisfy equation (8-11), 
This does not prove, of course, that R{p,po) = R{po,p). As it turns 
out, this is indeed equality, but since its validity can only be established 
by a rather lengthy analysis (Chandrasekhar), we will accept without 
proof 


( 8 - 12 ) 

as another manifestation of the principle of reciprocity. 

Equation (8-11) can be factored to give 

rf/j j^l + 2^io J d/i'j 

(8-13) 

Guided by the form of equation (7-87) we can define 

H{fi) ^l + 2n r R{n,ii') dll' (8-14) 

JQ 

and write equation (8-13) as 

and we see that equation (8-14) is another definition of the //-function. 
Equation (8-14) is exact in that it does not involve any orders of 
approximation. 

If we write equation (8-15) as 


(— + -) 
Vmo m/ 




4 (1 + IM) 


(8-16) 


and substitute equation (8-16) back into equation (8-14), we get 


= 1 + r (8-17) 

2 Jo p A- P^ 

which is the integral equation for H promised earlier. Equation (8-17) 
can be solved iteratively to determine H{p) to any degree of accuracy. 

Chandrasekhar presents a much more sophisticated derivation of 
equation (8-13) and the succeeding relations, resulting in an equation 
very similar to equation (8-11) for his scatter function; his equation 
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can, of course, be obtained from our equation (8-11) by using the 
correspondence equation (4-10). Once the physics of our derivation of 
equation (8-11) is fully understood, it can be of great benefit to review 
Chandrasekhar’s analysis, from the point of view of gaining facility in 
manipulating the fundamental definitions and using the integral form 
of the RTE to develop our results. 

The mean value of is a useful starting point for the 

iterative solution of equation (8-17). Define 


Ho= f dp 

Jo 

Multiply equation (8-17) by dp and integrate 


(8-18) 


Jo ^ Jo Jo M M 

Interchange p and /i' in the above and add the two results together 

"1 /•! H{p)H{p^)p' dp^ dp 


2/^'ifWd^ = 2 + |[/^7„ 

+/7 

JQ Jo 


fi + n' 

1 d/j, dfi'^ 

fi' + H 


= 2 + ^ f f dfi dfi' 

2 Uo Jo J 


or 


from which 


Ho = l + 


Ho = l[l- Vl-w) 


(8-19) 


We now write equation (4-17) for the planetary albedo for the model 
atmosphere in terms of the 7f- functions 


r(/xo) = 2 / R{p, po)p dp 

Jo 
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and using equation (8-16) 




f dp. 

^ Jo V MO + A*/ 

^H{po) f f H{p) dp - po r ^^dp 

^ Uo Jo + 


By use of equation (8-17) the integral can be eliminated to give, along 
with equation (8-19) 


r(/io) = 1 - H{po)y/l (8-20) 

Equation (8-20) is plotted in figure 8-2 for various values of a). The 
spherical albedo follows from equation (4-24) 

f - 1 - 2y/l - oj f pqH{p.q) dpo (8-21) 

Jo 

The first-moment integral in equation (8-21) was evaluated numerically 
from the H-function tables, and r vs. a) is plotted in figure 8-3. 

We can get an approximate analytic form for the /7-function, and 
thus the reflection coefficient from the two-stream solution. For the 
two-stream case, n = 1, p = l/v/3, a\ = 1, and the eigenvalues follow 
immediately from equation (7-59) 

k = y^3(l — Cj) 


From the definition of the //-function, equation (7-49) 


//(^) = 


1 “h fly/S 
1 “h — a}) 


and thus the reflection function becomes 


( 8 - 22 ) 


R{fi,po) 


^ (1 -h py/3){l + ^tQV^) 

4(/i + /io) [i + -(D)] [l + /iov^3(l - (D)j 


(8-23) 


If we evaluate H{pj for the n = 4 case used as our example, we get for 
fi = 0.5, Cj — 0.8, the data in table 8-1. 

It can be seen that going from a two-stream to an eight-stream model 
significantly improves the accuracy of evaluating the //-functions, and 
hence also the reflection functions. 
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TABLE 8-1. /io) 



2-Stream 

8-Stream 

Exact 

0.0 

0.53803 

0.56256 

0.56528 

0.1 

.49378 

.52749 

.53645 

0.2 

.46589 

.48915 

.49607 

0.3 

.44759 

.45399 

.45950 

0.4 

.43529 

.42289 

.42745 

0.5 

.42689 

.39559 

.39943 

0.6 

.42112 

.37152 

.37488 

0.7 

.41717 

.35022 

.35318 

0.8 

.41450 

.33124 

.33391 

0.9 

.41276 

.31423 

.31666 

1.0 

.41169 

.29891 

.30114 


First-Order Solution for the i/-Function 

The zeroth-order solution for H[ii) is given by equation (S-IO). If 
we use this in the right-hand side of equation (8-17), as the first guess 
of H{p)^ we get the first-order solution 


(m) = 1 + 'n (“^) (8-24a) 

A somewhat better approximation can be obtained by first solving 
equation (8-17) for H{p) to give 




1 _ g(/xQ y i ^ 
2 ^ io M . 


and then solving this by substituting i/o for H{p) on the right-hand 
side. This gives 




1 - In 



(8-24b) 


For Cj = 0.5, reflection coefficients computed from equation (8-20), using 
the first-order solutions for H from both equations (8-24a) and (8-24b), 
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are compared in figure 8-4 with the exact solution using the exact H- 
functions tabulated at the end of this chapter. Both approximations 
are adequate for large absorption (a; 1), but equation (8-24b) gives 

decidedly better results for larger to (nearly conservative scattering). 
Whether either approximation is adequate depends, of course, on the 
application. 

Equation (8-24a) or equation (8-24b) could, in turn, be resubstituted 
into equation (8-17) and a second-order solution derived. The algebra, 
however, becomes quite messy, and it is probably advisable to evaluate 
the resulting integrals numerically if this order of approximation is 
required. 

Table 8-2 compares results computed for equations (8-24a) and 
(8-24b) with the exact results. 

TABLE 8-2. VALUES OF i/i(/i) 



Equation (8-24a) 

Equation (8-24b) 

Exact 

0.0 

1.00000 

1.00000 

1.00000 

0.1 

1.07935 

1.07554 

1.07241 

0.2 

1.11315 

1.11727 

1.11349 

0.3 

1.13281 

1.14790 

1.14391 

0.4 

1.14615 

1.17202 

1.16800 

0.5 

1.15721 

1.19174 

1.18776 

0.6 

1.16848 

1.20826 

1.20436 

0.7 

1.18146 

1.22237 

1.21858 

0.8 

1.19706 

1.23459 

1.23091 

0.9 

1.21578 

1.24528 

1.24171 

1.0 

1.23785 

1.25473 

1.25128 


The iteration procedure for computing H{p) uses equa- 

tions (8-24), as the first guess and proceeds from equation (8-17). The 
solution converges fairly rapidly for small (D, but more and more iter- 
ations are needed as a) — ^ 1. Tables of the i7-function for isotropic 
scattering are included at the end of this chapter (table 8-3). 

If one uses equation (8-17) directly, the convergence proceeds some- 
what as shown in figure 8-5. 

Chandrasekhar, recognizing the slowness of the convergence of equa- 
tion (8-17), gives an alternate integral equation form for H[fj) 


1 

W) 



+ n' 


(8-25) 
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(8-24b) 
(8- 24a) 


Figure 8-4. Selected single-scatter eilbedo solution from figure 8-2 showing the 
accuracy of two of the approximate solutions, equations (8- 24a) and (8-24b). 
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H(y) 



Figure 8-5. Sketch illustrating the iterative behavior of equation (8-17). 



By a straight application of equation (8-25), however, no signif- 
icant improvement in the rate of convergence is noticed, although 
Chandrasekhar claims that it is decidedly superior to equation (8-17). 
Its convergence proceeds as sketched in figure 8-6. Convergence could 
perhaps be speeded up somewhat if we take, for example, the mean of 
the zeroth and first iterations as the second guess, the mean of the sec- 
ond and third iterations as the fourth guess, etc. This was not tried by 
the writer. Chandrasekhar’s iterative procedure is not discussed in his 
text, but perhaps this is the method he used to increase the rate of con- 
vergence. The problem is, of course, academic, as numerical solutions 
are available for all u), and the job is finished. 
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TABLE 8-3. iL-FUNCTIONS FOR ISOTROPIC SCATTERING 



w = 0.1 

1! 

p 

G) = 0.3 

U) = 0,4 

UJ = 0.5 

0.00 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

0.05 

1.007841 

1.016118 

1.024902 

1.034293 

1.044428 

0.10 

1.012385 

1.025632 

1.039895 

1.055387 

1.072402 

0.15 

1.015844 

1.032948 

1.051553 

1.071988 

1.094720 

0.20 

1.018645 

1.038919 

1.061149 

1.085784 

1.113465 

0.25 

1.020993 

1.043956 

1.069300 

1.097594 

1.129654 

0.30 

1.023006 

1.048296 

1.076365 

1.107899 

1.143889 

0.35 

1.024760 

1.052095 

1.082581 

1.117017 

1.156568 

0.40 

1.026306 

1.055459 

1.088110 

1.125169 

1.167971 

0.45 

1.027685 

1.058467 

1.093072 

1.132519 

1.178306 

0.50 

1.028922 

1.061177 

1.097559 

1.139192 

1.187734 

0.55 

1.030042 

1.063634 

1.101641 

1.145285 

1.196381 

0.60 

1.031060 

1.065875 

1.105375 

1.150876 

1.204347 

0.65 

1.031991 

1.067929 

1.108805 

1.156030 

1.211717 

0.70 

1.032846 

1.069820 

1.111971 

1.160799 

1.218559 

0.75 

1.033635 

1.071567 

1.114903 

1.165227 

1.224932 

0.80 

1.034365 

1.073187 

1.117626 

1.169352 

1.230885 

0.85 

1.035043 

1.074694 

1.120165 

1.173203 

1.236459 

0.90 

1.035674 

1.076099 

1.122536 

1.176810 

1.241693 

0.95 

1.036264 

1.077413 

1.124758 

1.180196 

1.246617 

1.00 

1.036816 

1.078645 

1.126844 

1.183380 

1.251259 

Mean 

1.026334 

1.055728 

1.088933 

1.127017 

1.171573 

1st mom. 

0.515611 

0.533155 

0.553122 

0.576214 

0.603486 
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TABLE 8-3. Continued 



II 

o 

II 

P 

II 

o 

00 

u) — 0,85 

(D = 0.90 

0.00 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

0.05 

1.055513 

1.067885 

1.082180 

1.090455 

1.099980 

0.10 

1.091388 

1.113078 

1.138860 

1.154176 

1.172201 

0.15 

1.120454 

1.150357 

1.186654 

1.208633 

1.234933 

0.20 

1.145168 

1.182519 

1.228642 

1.257015 

1.291436 

0.25 

1.166734 

1.210934 

1.266321 

1.300861 

1.343268 

0.30 

1.185867 

1.236418 

1.300586 

1.341086 

1.391346 

0.35 

1.203043 

1.259515 

1.332031 

1.378299 

1.436276 

0.40 

1.218599 

1.280617 

1.361086 

1.412937 

1.478491 

0.45 

1.232788 

1.300016 

1.388077 

1.445334 

1.518322 

0.50 

1.245806 

1.317943 

1.413259 

1.475753 

1.556029 

0.55 

1.257807 

1.334580 

1.436839 

1.504405 

1.591821 

0.60 

1.268919 

1.350077 

1.458986 

1.531467 

1.625876 

0.65 

1.279244 

1.364560 

1.479845 

1.557089 

1.658340 


1.288870 

1.378134 

1.499537 

1.581397 

1.689343 

0.75 

1.297870 

1.390887 

1.518166 

1.604501 

1.718996 


1.306306 

1.402899 

1.535825 

1.626499 

1.747398 

0.85 

1.314234 

1.414235 

1.552593 

1.647475 

1.774634 


1.321700 

1.424955 

1.568540 

1.667505 

1.800784 

0.95 

1.328745 

1.435109 

1.583730 

1.686656 

1.825916 


1.335406 

1.444745 

1.598217 

1.704989 

1.850095 

Mean 

1.225148 

1.292221 

1.381966 

1.441651 

1.519494 

1st mom. 

0.636634 

0.678670 

0.735817 

0.774378 

0.825317 
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TABLE 8-3. Continued 



Cj = 0.92 

(D = 0.94 

(D = 0.96 

u? = 0.98 

a) = 0.99 

0.00 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

0.05 

1.104330 

1.109161 

1.114731 

1.121700 

1.126408 

0.10 

1.180588 

1.190025 

1.201077 

1.215196 

1.224940 

0.15 

1.247339 

1.261434 

1.278137 

1.299801 

1.314988 

0.20 

1.307860 

1.326672 

1.349186 

1.378764 

1.399768 

0.25 

1.363707 

1.387286 

1.415750 

1.453571 

1.480738 

0.30 

1.415788 

1.444169 

1.478699 

1.525056 

1.558707 

0.35 

1.464702 

1.497907 

1.538598 

1.593750 

1.634178 

0.40 

1.510876 

1.548914 

1.595842 

1.660019 

1.707496 

0.45 

1.554634 

1.597502 

1.650726 

1.724130 

1.778906 

0.50 

1.596228 

1.643917 

1.703480 

1.786286 

1.848594 

0.55 

1.635867 

1.688357 

1.754287 

1.846650 

1.916703 

0.60 

1.673720 

1.730986 

1.803301 

1.905354 

1.983349 

0.65 

1.709935 

1.771944 

1.850652 

1.962507 

2.048627 

0.70 

1.744637 

1.811351 

1.896449 

2.018205 

2.112617 

0.75 

1.777935 

1.849314 

1.940791 

2.072528 

2.175388 

0.80 

1.809926 

1.885924 

1.983763 

2.125549 

2.236997 

0.85 

1.840696 

1.921266 

2.025441 

2.177330 

2.297498 

0.90 

1.870323 

1.955413 

2.065896 

2.227929 

2.356936 

0.95 

1.898876 

1.988434 

2.105188 

2.277398 

2.415353 

1.00 

1.926417 

2.020389 

2.143376 

2.325784 

2.472787 

Mean 

1.559038 

1.606492 

1.666667 

1.752201 

1.818182 

1st mom. 

0.851467 

0.883087 

0.923548 

0.981749 

1.027182 
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TABLE 8-3. Concluded 



Cj = 0.995 

Cj = 0.999 

0.00 

1.000000 

1.000000 

0.05 

1.129618 

1.133736 

0.10 

1.231690 

1.240491 

0.15 

1.325628 

1.339664 

0.20 

1.414627 

1.434417 

0.25 

1.500122 

1.526160 

0.30 

1.582903 

1.615664 

0.35 

1.663460 

1.703400 

0.40 

1.742120 

1.789679 

0.45 

1.819117 

1.874719 

0.50 

1.894621 

1.958680 

0.55 

1.968765 

2.041679 

0.60 

2.041654 

. 2.123810 

0.65 

2.113372 

2.205145 

0.70 

2.183989 

2.285743 

0.75 

2.253563 

2.365652 

0.80 

2.322145 

2.444913 

0.85 

2.389778 

2.523558 

0.90 

2.456500 

2.601616 

0.95 

2.522344 

2.679111 

1.00 

2.587341 

2.756066 

Mean 

1.867918 

1.938693 

1st mom. 

1.061731 

1.111331 


156 



Chapter 9 


Additional Topics 

There are a great number of problems contained in radiative transfer 
theory that were not addressed at all in these notes. We will mention 
just a few of these here as a conclusion to the text, describe them 
briefly, and indicate some references which perhaps address them more 
thoroughly. 

Determination of Optical Parameters 

All of the methods discussed in the text have assumed that the 
optical parameters used in the equations, such as optical depth, phase 
function, asymmetry parameter, single-scattering albedo, etc., were 
all known. These parameters can be computed to an acceptable 
degree of accuracy, in most cases, by the use of well established 
numerical or theoretical methods or both. The complete repertoire of 
procedures again consists of both “exact” and approximate methods, 
but unfortunately, it would take another text larger than the present 
one to describe them in sufficient detail. 

For homogeneous atmospheres, the optical depth can be computed 
for a single frequency and along a given slant path with little difficulty. 
Unfortunately, all measuring devices measure radiation in a finite band 
of frequencies, with a variable response across the band. The absorption 
coefficient varies very rapidly with frequency, and at a single frequency, 
many tens of nearby lines may contribute to the monochromatic ab- 
sorption coefficient. Thus, a great deal of data concerning the positions 
of line centers, line strengths, and line shapes must be available. Addi- 
tionally, since the lines are generally spaced such that they overlap to 
varying degrees, a very fine grid spacing in wavelength must be used 
to get the total absorption in a given finite bandwidth. This prob- 
lem is further complicated by the fact that the line optical parameters 
vary strongly with altitude (i.e., pressure and temperature) and with 
frequency, and hence, the absorption changes in a strongly nonlinear 
fashion with these parameters. As a result, generally three nontrivial 
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integrations are required to completely describe the absorption charac- 
teristics of radiation; i.e., over wavelength, angle, and altitude. 

Band models attempt to reduce the frequency integration to a 
tractable problem. Some assumptions concerning the distribution of 
line centers and the distribution of line strengths are made to reduce the 
frequency integral to one which can be evaluated in terms of elementary 
functions. This scheme has produced a number of popular band models 
which have been used in a number of atmospheric physics applications, 
such as climate modeling and studies of the thermal structure of the 
atmosphere. Three excellent references for band model derivations and 
applications are those by Goody (1964), Rodgers (1976), and Anding 
(1969). 

The scattering optical properties can also be determined, at least for 
spherical particles. The scatter properties of particles which are very 
small relative to the wavelength of the incident radiation (i.e., molecular 
scattering) can be accurately described by the Rayleigh theory (Liou, 
1980, van de Hulst, 1957), while for very large particles, ray tracing 
techniques are generally used (Liou, 1980). Particles in the intermediate 
size range are the ones which cause most of the computational problems. 
The most general theory here is the Mie theory (Liou, 1980, van de 
Hulst, 1957, Stratton, 1947). However, as stated earlier, this theory 
is complete only for spherical particles, although some success with 
cylinders and flat plate particles has been reported (Liou, 1980, van 
de Hulst, 1957, Kerker, 1969). The optical properties such as phase 
function, and scatter and absorption cross sections are functions of only 
two parameters — the ratio of the wavelength of the incident light to the 
particle radius, and the complex index of refraction of the material 
from which the particles are made — and are independent of pressure 
and temperature. 

If the particle size distribution is known (number density of particles 
as a function of particle radius, for example), the overall properties of a 
unit volume of scatterers (polydispersion) can be computed (Liou, 1980, 
Deirmendjian, 1969). Since both the total number of particles per unit 
volume and their size distribution may in general vary with altitude, 
there is thus a strong altitude dependence built into the scattering 
properties of a polydispersed conglomerate of particles. 

As mentioned above, these computational procedures are well de- 
fined and are used extensively in the literature, although the computa- 
tional details are quite involved and time consuming, and in some cases 
tax even the most modern of high-speed computers. 

Some approximations to the Mie results have been reported in the 
literature, and may be profitably used in studies in which the ultimate 
in accuracy is not needed — e.g., in studies of climate modeling and the 
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effects of aerosols on global climate (van de Hulst, 1957, Penndorf, 1962, 
Plass, 1966). 

Finite Homogeneous Atmospheres 

The discrete ordinates method (chap, 7) and the principle of invari- 
ance (chap. 8), as discussed in the text, are applicable only to semi- 
infinite homogeneous atmospheres with isotropic scattering. Chan- 
drasekhar (1960) and Sobolev (1975) extend these techniques in elegant 
mathematical fashion to finite homogeneous atmospheres with arbitrary 
(to some extent) phase functions. 

The principle of invariance can be simply stated for a finite atmo- 
sphere. If we add an infinitely thin layer of the same optical properties 
to the top of a finite atmosphere, then the changes in the reflection func- 
tion and the transmission function for this incremented atmosphere can 
be computed. Similarly, if we add an infinitely thin layer to the bot- 
tom of the atmosphere, these changes can again be computed. The 
two sets of change must be equal; equating them, one arrives at two 
coupled nonlinear integro-differential equations for the reflection and 
transmission functions for the finite atmosphere, equations quite simi- 
lar to equation (8-10). By rearranging terms and factoring as we did 
earlier, two functions can be defined which are similar in utility to 
the //-functions. These are the famous X- and F-functions of Chan- 
drasekhar, which describe the reflection and transmission of isotropic 
radiation in finite atmospheres. Further manipulation yields a coupled 
set of integral equations for these functions similar to equation (8-17). 
Then the angular distribution of the radiant energies from both the top 
and the bottom of the atmosphere can be described in equations similar 
to equation (7-84). 

The X- and F-functions can be computed and tabulated (see 
Chandrasekhar), just as we did for the //-functions, for specific phase 
functions. 

As the thickness of the finite atmosphere increases, approaching the 
semi-infinite case, the X-function approaches the //-function and the 
F-function goes to zero. Thus, the X-function is related to the reflection 
properties of the finite atmosphere, while the F-function is related to 
its transmission properties. 

For very thin atmospheres, the X-function approaches unity and 
the F-function approaches and the resultant equations for the 

radiance reduce to the single-scattering solution we found in chapter 5. 

Anisotropic Scattering 

In chapter 1 it was pointed out that the phase function can in many 
cases be described by a Legendre polynomial expansion in the scattering 
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angle. In both the semi-infinite and finite atmospheres, if this is done 
and the principle of invariance applied, the result is one (for semi-infinite 
atmospheres) or two (for finite atmospheres) integral equations for an 
/f-function or for the X- and F-functions, for each term of the Legen- 
dre expansion, and these are in general horribly coupled. The numerical 
problems thus generated are so enormous, that, to this writer’s knowl- 
edge, no one has generated general tables of these functions except for 
isotropic scattering and some limited results for Rayleigh scattering. 
However, even for the relatively simple cases of the two-term expansion 
and the two- term Rayleigh expansion, the computational difficulties 
are such that even Chandrasekhar only presents a limited number of 
numerical tables. 

Similarity may be applied in some cases. However, in general some 
other numerical technique, such as the adding or doubling methods 
to be described later, is usually used. The discrete ordinates method, 
and the related spherical harmonics methods, have been successfully 
applied numerically in some limited cases, even for nonhomogeneous 
atmospheres, but it is generally conceded that the other procedures are 
numerically and computationally superior for these applications. 

Effect of Surface Albedo 

In our analysis of the inclusion of surface effects in chapter 4, it 
was assumed that the reflective surface was Lambertian (i.e., isotropic 
scatter from the surface) and was the same for all parts of the surface 
plane. This is in general not a realistic approximation, but again the 
numerical results may be accurate enough for some applications. Little 
work has been done on other than Lambertian surfaces, but some results 
are available for specularly reflecting surfaces. See the thesis by Tanre 
and the paper by Deschamps et al. in Deepak, 1980. 

Other Computational Techniques 

There are a number of so-called “exact” methods available in the 
literature, which are comparable to or somewhat better than the 
discrete ordinates method for nonisotropic scattering. These methods 
are exact in the sense of some limiting process as described with the 
discrete ordinates method covered in chapter 7. A few of these methods 
will be discussed in this subsection. 

Adding and doubling methods. These methods are similar and are 
both based on the following premise: suppose we have two slabs of opti- 
cally active material, and suppose that we know the reflective and trans- 
mission properties of each slab separately. If we place the two slabs to- 
gether, face to face, then by considering the multiple transmissions and 
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reflections between the two slabs, it is possible to determine the overall 
transmission and reflective properties of the composite slab considered 
as a unit. In its most fundamental form, we can show this as follows 
in figures 9-1 and 9-2, where the slabs are shown separated for clarity 
only. We let the reflectance and transmission coefficients be denoted by 
R and T, respectively, with subscript 1 referring to the upper slab and 
subscript 2 referring to the lower. 



Figure 9-1. Illustration of the various orders of scattering between two finite thick 
layers. The layers are shown separated for clarity only. 


The following rays emerge from the top of the composite slab: 

1. Ray 1 is simply the reflection from the upper slab, R\. 

2. Ray 2 is a ray transmitted through the upper slab, reflected from 
the lower slab, and transmitted through the upper slab, T 1 R 2 T 1 . 

3. Ray 3 is transmitted through 1, reflected from 2, reflected back 
down from 1, reflected again from 2, and transmitted out through 

1, TiR2^1^2^1- 

4. For the remaining rays, there are similar multiple reflections between 
1 and 2 and transmissions through 1. 

Collect all these together, and we have for the total reflection from 
the top of layer 1 

R\2 — H- T 1 R 2 T 1 + + Tii?2^1^2^1^2^1 -h ’ • • 

and these can be collected to give 

Ri 2 = Ri+ TiR2Ti{l + R 1 R 2 + RiH^ + ■ ■ •) 

Since R 1 R 2 < i we can write this last as 


R\2 = R\ 


RlTf 
I — RIR2 


(9-1) 
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Transmission is handled the same way, as shown in figure 9-2, and the 
composite transmission function can be written 

T \2 = T1T2 -l- T1R2R1T2 + T1R2R1R2R1T2 
^TiT 2 [l + R1R2 + RjRi + ---) 



Figure 9-2. Same as figure 9-1 except showing the various orders of scattering 
involved in diffuse transmission. 


or 


Ti2 


T1T2 
1 — R1R2 


(9-2) 


The similarity between equation (9-1) and equation (4-39) or equa- 
tion (4-40) cannot have escaped the reader’s attention. 

In more realistic application, the order of multiplication in equa- 
tions (9-1) and (9-2) must be preserved, and the simple products are 
replaced by integral functions over all directions (see, e.g., Liou, 1980). 
Another approach is to construct the R and T as matrices whose ele- 
ments are in general integrals of sundry combinations of the directional 
representations of the reflection and transmission coefficients. This ap- 
proach is directly oriented toward computer application; see, for exam- 
ple, the excellent paper by Twomey, Jacobowitz, and Howell (1966); 
see also van de Hulst (1963), the highly mathematical series of papers 
by Grant and Hunt (1968a, 1968b, 1969), and the paper by Hunt and 
Grant (1969). 

Both the adding and the doubling methods use the generalized form 
of equations (9-1) and (9-2). 

The doubling method is applicable to homogeneous atmospheres. 
A very thin slab is selected, say Ar = 2 x 10“^^ or so, and the 
reflection and transmission coefficients are computed by one of the thin- 
atmosphere solutions covered in chapter 5, say the thin-atmosphere 
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solutions, equations (5-16) and (5-19). Now, if we assume two slabs 
of the same thickness and same R and T, then the generalized forms 
of equations (9-1) and (9-2) can be used to compute R and T for the 
thickness 2 At. These R and T then can be used in turn to compute a 
pair of R and T for a thickness 4 At, 8 At, etc., doubling the thickness 
with each application of equations (9-1) and (9-2) until the desired total 
optical thickness is reached. Obviously, the doubling method can only 
be used for homogeneous optical slabs. 

The adding method is similar and can be used for inhomogeneous 
atmospheres. Suppose the inhomogeneous layer is divided into a 
number of thin layers, and the thin-atmosphere solution (or indeed, 
the doubling method!) is used to compute the R and T for each layer 
separately. Then, the appropriate forms of equations (9-1) and (9-2) 
can be used to get R \2 and T\ 2 - Then the third layer can be added to 
this to yield R\2^ and T123, and then the fourth and succeeding layers 
until the entire atmosphere is completed. Liou (1980) presents some 
very useful tables of reflection and transmission coefficients computed 
from both the discrete ordinates method and the doubling method. 

Coakley, Cess, and Yurevich (1983) present an interesting method of 
using reflection and transmission coefficients computed from the delta- 
Eddington method with the adding and doubling methods. 


The spherical harmonics method. The spherical harmonics method 
is very similar to the discrete ordinates method; in fact, the discrete 
ordinates method is a specialized form of the spherical harmonics 
method. In this method, shown here for the azimuthally symmetric 
case, it is assumed that the intensity itself as well as the phase function 
can be expanded in a series of Legendre polynomials 




V— ' 2m T 1 n \ / r \ 

m=0 


(9-3) 


where the ^m(^) are coefficients which are functions of r only, and the 
phase function is expanded as 


N 

= ^{2n+l)fnPn{f^)Pn{fi’) 


n=0 


(9-4) 


163 



Introduction to the Theory of Atmospheric Radiative Transfer 

After substituting into the RTE and simplifying, one gets a system 
of differential equations for the 

(m+ 1) — + (2m - 1)(1 -w/m)0m 

= 47t(1 - w)B(T)<5o”^ (/o = 1) (9-5) 

where B{T) is the Planck function. For isotropic scattering, all the 
fm — 0 except for /g = 1. If one retains N terms of the expansion 
equation (9-3) (the Pn-approximation), then equation (9-5) yields a set 
of iV -|- 1 simultaneous, linear ordinary differential equations for the 
t/’O) V’l^ •j (V'n+l is set to zero), which when substituted back 
into equation (9-3), give the intensity at any r,/i. 

The boundary conditions for the P„-method are difficult to satisfy 
exactly, and in general some approximation must be used. We have seen 
this earlier in connection with the Eddington method, and in fact the 
spherical harmonics method with iV = 0 and iV = 1 yields precisely the 
Eddington equations. See Ozisik (1973) for an elementary discussion 
of the P„-method and for some schemes for handling the boundary 
conditions. See also the discussions in Kourganoff (1963) and Lenoble 
(1977). 

Monte Carlo method. This is perhaps the only method known 
which can be applied to any radiative transfer problem regardless of 
asymmetry, nonhomogeneity, or any other anomaly, and is the only 
method which can really be called “exact.” However, as in other 
applications, there is no “free lunch,” and one must pay a heavy price 
in computer costs — mostly time — for this flexibility and general utility. 

Basically, in the Monte Carlo method one injects a series of single 
photons into the medium and follows one photon at a time in space and 
time as it travels through the three-dimensional medium. Whenever 
the photon encounters an absorber or scatterer, a suitable probability 
is used to determine whether an actual interaction occurs and what 
type. If the interaction is an absorption, the computations stop here 
and the energy of the photon is used to increment the total energy of the 
medium, and another photon is injected into the medium and followed. 
If the interaction is a scattering, the direction into which the photon is 
scattered is determined probabilistically from the phase function. The 
photon is followed through ensuing scatterings or absorptions or until 
it escapes through the top of the atmosphere. It is apparent that a 
great many photons must be tracked (orders of hundreds of thousands) 
to provide a reliable sample size from which to determine reflection, 
transmission, and absorption distributions, and therefore a great deal 
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of computer time is required. A number of computational schemes have 
evolved to shorten the computational time and retain the accuracy of 
this method, but the expense has precluded its wide application for 
radiative transfer studies. Its utility seems to be in the areas where 
absolutely nothing else works, and to provide some limited benchmark 
results against which to compare the results of more rapid but perhaps 
less precise analyses. 

See the most interesting discussions of this method in Irvine and 
Lenoble (1973) and in the paper by Hansen and Travis (1974). The 
papers by Kattawar and Plass (1968) and Plass and Kattawar (1968) 
best describe the application of this method to radiative transfer 
problems. 

There are, of course, many other methods not mentioned here 
for solving the RTE to greater or lesser degrees of approximation. 
These include the method of successive orders of scattering, which is 
an extension of the single-scatter method derived in chapter 5, the 
eigenvalue expansion method of Case, the Gauss-Seidel method (a 
numerical technique), and many others. Some of these are briefly 
discussed in Irvine and Lenoble (1973), where specific references are 
given, as well as in the paper by Hansen and Travis (1974), and the 
text by Ozisik (1973). A much more comprehensive discussion is given 
in Lenoble (1977), with many references and the basic equations. 

Non-Homogeneous Atmospheres 

Practically all of the methods discussed in the present text are 
restricted to solutions in a homogeneous atmosphere. There has been 
much more effort expanded in applying these methods approximately 
to nonhomogeneous atmospheres. What is generally done is to divide 
the atmosphere into a number of thin layers and treat each layer by, for 
example, the discrete ordinates method. The main difference between 
this technique and the approach we have taken is in the application of 
the boundary conditions. Here, one can use the zero diffuse radiation 
boundary condition only at the top of the uppermost layer and at the 
bottom of the lowest layer. In between, the boundary conditions must 
be set up to insure continuity of flux or energy across each boundary. 
This procedure generally leads to a set of simultaneous algebraic 
equations which must be solved for a set of constant coefficients for 
each layer considered (e.g., the constants A and B of equations (5-69) 
and (5-70) would be different in each layer). Liou (1973) has done this 
using the discrete ordinates method, and IViscombe (1977) has used 
the delta-Eddington method in the same way. Comparisons with the 
more nearly exact adding method indicate that good accuracy can be 
obtained with these schemes. 
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Other Problems 

Finally, there are many other problems in radiative transfer theory 
which are seldom mentioned in the literature. There is, for example, 
the inclusion of horizontal inhomogeneities in the plane-parallel atmo- 
sphere we have been using (e.g., finite clouds or actual differences in 
optical properties due to climatic or meteorological effects), or the re- 
vision of the plane-parallel assumption itself, i.e., considering spherical 
atmospheres. Some of these problems have been addressed by neutron 
physicists, since the travel of neutrons through absorbing and scattering 
media is described by an equation very similar to our radiative transfer 
equations— the main difference being that the neutrons can travel with 
different speeds, while our photons all travel at the speed of light. 

Other problems include the shadow effect, in which the shadowing 
of one particle by another prevents the second particle from interacting 
fully with the incident field— it is shielded to some extent by the 
particles in front of it. This can occur, according to van de Hulst 
(1957), if the mean spacing between particles becomes less than four 
or five particle diameters. This problem practically never arises in 
atmospheric applications of radiative transfer theory, but can arise in 
neutron theory. 

Another major problem area, which just over the last ten years or 
so has begun to receive attention in the literature, is the problem 
of the transfer of polarized radiation components and their use in 
studying the properties of atmospheric components, and particularly 
in the study of radiation from the surface of the oceans and clouds. 
In most cases, this can be handled both numerically and analytically 
by replacing the scalar equation we have been using with a vector 
equation; i.e., the intensity scalar becomes a four-component vector 
whose components are usually the parameters (see Deirmendjian, 

1969, van de Hulst, 1957, Hansen and Travis, 1974), and the phase 
function becomes a 4 x 4 phase matrix, whose components characterize 
the polarization produced by a single act of scattering. Many of the 
numerical techniques discussed in this chapter (adding, doubling, etc.) 
can be used to analyze polarized fields, but little analytic work has 
been done in this area (see Irvine and Lenoble, 1973, and Lenoble, 
1977). Polarization for Rayleigh scattering has been considered by 
Chandrasekhar (1960), and some work has been done by Sekera (see 
Lenoble, 1977). 
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